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Many superconducting materials are composed of weakly coupled conducting layers. Such a 
layered structure has a very strong influence on the properties of vortex matter in a magnetic field. 
This review focuses on the properties of the Josephson vortex lattice generated by the magnetic field 
applied in the layers direction. The theoretical description is based on the Lawrence-Doniach model 
in the London limit which takes into account only the phase degree of freedom of the superconducting 
order parameter. In spite of its simplicity, this model leads to an amazingly rich set of phenomena. 
We review in details the structure of an isolated vortex line as well as various properties of the vortex 
lattice, both in dilute and dense limits. In particular, we present an extensive consideration on the 
infiuence of the layered structure and thermal fluctuations on the selection of lattice configurations 
at different magnetic fields. 



uniaxial continuous superconductor breaks down, which 
happens when the layer separation d is greater than the 
z-axis superconducting coherence length, d ^ ^c- 

The most prominent example is the high-Tc cuprate su- 
perconductors, discovered in 1986^"—, which led to a huge 
interest in physics of layered superconductors. The two 
most studied cuprate compounds, YBa2Cu307 (YBCO) 
and Bi2Sr2CaCu20a; (BSCCO), have similar transition 
temperatures T^ ~ 90K and represent two important 
particular cases. YBCO is moderately anisotropic, with 
the anisotropy factor 7 sa 5 — 7, and its "layeredness" 
becomes essential at low temperatures when the c-axis 
coherence length ^c drops below the layer spacing d. On 
the other hand, BSCCO has a huge anisotropy factor, 
7 « 400 — 1000, and behaves as a layered supercon- 
ductor practically in the whole temperature range below 
Tc- Other naturally layered superconductors include the 
transition metal dichalcogenides^'^ and organic charge- 
transfer salts formed with the molecule BEDT-TTF— i^. 
Important new family of atomically layered supercon- 
ducting materials, iron pnictides and chalcogenides, has 
been discovered in 2008^- and are being extensively ex- 
plored since then, see, e.g., Reviewaii"— . Anisotropy 
of most compounds is actually not very high and typi- 
cally they behave as anisotropic three-dimensional ma- 
terials. There are, however, important exceptions. The 
most studied compound in which the layered structure 
is clearly essential is SmFeAsOi_a;FaJ^ with T^ up to 
bbK. For example, Josephson nature of the in-plane 
vortices at low temperatures has been recently demon- 
strated in this compoundi^. Also, several iron pnictide 
compounds with extremely high anisotropy have been 
discoveredi^^— . Properties of these compounds remain 
mostly unexplored due to their rather complicated com- 
position. 

All layered superconductors share a very similar gen- 
eral behavior of the vortex matter generated by external 



FIG. 1. Illustration of a dilute lattice of Josephson vortices 
generated in a layered superconductor by magnetic fleld ap- 
plied along the layer direction. 



I. INTRODUCTION 

Layered superconductors are materials made from a 
stack of alternating thin superconducting layers sepa- 
rated by non-superconducting regions. The supercon- 
ducting layers are essentially two-dimensional (2D) as 
long as they are so thin that there is no variation in 
fields, or in the superconducting order parameter, across 
each layer. Such structures frequently occur naturally 
in anisotropic crystals. A layered superconductor can 
carry supercurrents along the layers, as well as between 
the layers. This is due to the Josephson tunneling of 
Cooper pairsi across the insulating regions that sepa- 
rate neighboring superconducting layers, i.e., each pair 
of neighboring layers forms one Josephson junction. In 
general, the z-axis (Josephson) supercurrents are weaker 
than the supercurrents along the layers. A mere "lay- 
eredness" of atomic structure, however, does not make a 
material automatically a layered superconductor. When 
the interlayer electrical coupling is strong enough, this 
discrete system of layers approximates to a continuous 
superconductor with uniaxial anisotropy. The case we 
are interested in here is when the approximation to a 



magnetic field, which is insensitive to the microscopic 
nature of superconductivity inside the layers. Several ex- 
cellent review articles have been published in the past 
covering different aspects of the vortex matter in type-II 
superconductorsi^r— . Nevertheless, we feel that further 
progress on the understanding of the Josephson vortices 
in layered superconductors warrants a specialized review, 
which will provide much more details and discuss impor- 
tant results obtained in recent years. 

This short review narrowly focuses on the vortex lat- 
tice which appears at magnetic fields applied along the 
layers. In this case the flux line winds its phase around 
an area between two neighboring layers and is called a 
Josephson vortex in analogy with a vortex in a supercon- 
ducting tunneling junction. The Josephson vortex con- 
tains out-of-plane currents that tunnel via the Joseph- 
son effect from layer to layer. The current distribution 
around vortex is anisotropic. As a consequence, the vor- 
tex lattice is also anisotropic, it is a triangular lattice 
strongly stretched along the layers, see Fig. [T] In ad- 
dition, the restriction to lie between the layers leads to 
commensurability effects and an energy barrier to tilt- 
ing the field away from the layers. There are two very 
different regimes depending on the strength of the mag- 
netic field Bx- The crossover field scale, i?cr, separat- 
ing these two regimes is set by the anisotropy factor 7 
and the layer periodicity, d, Bcr = $o/(27r7(i^), where 
<f>o = hc/2e is the flux quantum. In the case of BSCCO 
this field scale is around 0.5 tesla. In the dilute lattice 
regime, B^ < Bcr, the nonlinear cores of Josephson vor- 
tices are well separated and the distribution of currents 
and fields is very similar to that in continuous anisotropic 
superconductors^^. The dense lattice regime is realized 
at high fields, B^ > Ba, where the cores of the Josephson 
vortices overlap. In this regime the Josephson vortices fill 
all layers homogeneously^'^. This state is characterized by 
rapid oscillations of the Josephson current and by very 
weak modulation of the in-plane current. In this review 
we characterize in more detail these two lattice regimes. 

Note that we do not consider in this review properties 
of vortices generated by magnetic field applied perpen- 
dicular to the layers, along the c axis.^^ The structure 
of a c-axis vortex is very different from the structure of 
an in-plane vortex. In layered superconductors the c-axis 
vortex can be viewed as a stack of weakly coupled point- 
like pancake vortices. Properties of the pancake vortex 
lattice were also extensively explored, see, e.g.. Reviews 
[23I and [23 and references therein. 

Several experimental techniques have been employed 
to explore the Josephson vortex lattices. The dilute 
stretched lattice at small fields (< lOOG) has been di- 
rectly observed in YBCO with Bitter decoration by Dolan 
et al^, where the elliptical distribution of the flux 
around each Josephson vortex was also seen. At high 
flelds (> 1 tesla) the commensurability between the c- 
axis parameter of the Josephson vortex lattice and the 
interlayer separation leads to magnetic held oscillations 
which have been observed experimentally in underdoped 



YBCO in irreversible magnetization22i2£ and nonlinear 
resistivity"^^. 

In much more anisotropic BSCCO direct observation 
of Josephson vortices is not possible. However, when 
the magnetic held tilted at small angles with respect 
to the layers, the c-axis held component generates the 
pancake-vortex stacks which preferably enter the super- 
conductor along the Josephson vortices forming chains. 
Visualizing the flux of these chains, it is possible to find 
locations of vertical rows of the Josephson vortices and 
measure the in-plane lattice parameter ay. This was 
done using a variety of visualization techniques, such as 
Bitter decoration o^^'^'^ , scanning Hall probes^, Lorentz 
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microscop y^ '-'i'^" and magnetooptical imagingaiiHS These 
observations have been summarized in the Review |39J. 

Most extensively, properties of the Josephson vortex 
lattice have been explored in BSCCO using c-axis trans- 
port in small-size mesas^"— . These studies revealed a 
very rich dynamical behavior of the lattice, which is be- 
yond the scope of this review. The very important feature 
is that, due to low dissipation, the Josephson vortex lat- 
tice can be accelerated up to very high velocities. It is 
clear that understanding dynamics is not possible with- 
out good understanding of static lattice properties. The 
dynamic phenomenon closely related to static lattice con- 
figurations is magnetic-field oscillations of resistance for 
very slow lattice motion which have been discovered and 
explored in small-size BSCCO mesas^"— . The oscilla- 
tion period may correspond to either flux quantum or 
half flux quantum per junction depending on the mag- 
netic held and lateral size of the mesa. Interplay be- 
tween the bulk shearing interaction and the interaction 
with edges leads to very nontrivial evolution of lattice 
structures which we consider in this review. 

This review is organized as follows. We start in Sec. 
im in which we present the energy functional and equi- 
librium equations for the phase and vector potential. In 
Sec, mil we describe the structure and energetics of a sin- 
gle flux line. In Sec. |IV] we discuss the dilute JVL and 
consider in details the role of layered structure in select- 
ing lattice configurations. The properties of the dense 
JVL at high fields are considered in Sec.|Vl In this regime 
the structure and energy of the lattice can be evaluated 
analytically using expansion with respect to the Joseph- 
son coupling. In this section we also review the magnetic 
field dependence of lattice configurations and oscillations 
of the critical current in finite-size samples. Elastic prop- 
erties of both dilute and dense lattice are discussed in the 
corresponding sections. Based on the elastic energies, in 
Sec. |Vl]we review effects caused by thermal fluctuations. 



II. ENERGY FUNCTIONAL AND EQUATIONS 

FOR THE SUPERCONDUCTING PHASES AND 

VECTOR-POTENTIAL 

Theoretical analysis of the Josephson vortex matter in 
layered superconductors is based on a phenomenologi- 



cal model in which only phase degree of freedom of the 
superconducting order parameters is taken into account 
and its amplitude variations are neglected. 
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+ -pi^- C0s((/)„+i - 0„ + Xn.n+l)] I , (1) 

where Eq = ^Qd/{16TT^X^^^) gives the in-plane phase stiff- 
ness and £■,] = Eq/j'^ = <i>gd/(167r^A^) is the phase stiff- 
ness for smooth inter-layer phase variations, Aab and Ac 
are the components of the London penetration depth, 
and 7 = Ac /Aab is the anisotropy factor. The z- 
componcnt of the vector potential enters the tunneling 

term in the form Xn,n.+i = (2e//ic) J^^ dz Az^ Near 
the transition temperature the above phase model can be 
obtained from the celebrated Lawrence-Doniach model^ 
by fixing the order-parameter amplitude (London ap- 
proximation). However, the model is actually more gen- 
eral and describes Josephson properties of layered mate- 
rial in the whole temperature range. Starting from the 
phase model, a rich variety of the lattice properties can 
be derived, which we review in this article. 

Subject to some given boundary conditions, the config- 
uration of {(pm A} is determined by minimizing the free 
energy. This leads to a set of differential equations, e.g. 
minimizing with respect to the phase gives the current- 
conservation condition. 
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with the gauge-invariant phase difference defined as 
^n,n+i=4'n+i-4>n + Xn,n+i- In this cquatiou the Joseph- 
son length, Aj = ^d, appears for the first time. This 
length plays a very important role in layered supercon- 
ductors, as it determines the scale over which the phase 
can relax to minimize the Josephson coupling energy 
without costing too much energy in the gradient term. 
Three more equations result from minimizing with re- 
spect to the three components of the vector potential. 
We can write these in terms of the electric current density 
by using the Maxwell equation j ~ (c/47r)V x (V x A), 
giving. 
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where Jj|_„ is the 2D current density in the nth. layer and 
Jn^n+i is the current density in the z-direction between 
the nth and (n-l-l)th layer which has the maximum value. 
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The four equations (l2])-(|4|), are the starting point for 
finding the structure of vortices in layered superconduc- 
tors. In fact, we can make the job of solving this set of 
equations slightly clearer by combining them to find a 
differential equation for the gauge-invariant phase differ- 
ences alone. This is done by using the general result. 

And rin+i)d 

Jn,n+1 = / dz [V X (V X A)] 3 

C Jnd 

= M • (A||,„+i--A||^„)-— \7| Xn,n+1, (6) 

and combining this with ([2]) and (|4]) to arrive at, 

^|^Vn,n+l + ^2 ^^^ ^n.n+l (7) 
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[sin(^„+i,„+2-2sin(/j„,„+i+sin(^„_i^„] = 0. 
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Solving this equation will give the entire solution for cur- 
rents using (Ul to find Jn.n+ii and the conservation law. 
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III. STRUCTURE OF A JOSEPHSON VORTEX 
IN A LAYERED SUPERCONDUCTOR 

If we place a flux line directed along the layers, the sin- 
gularity associated with the vortex core can be avoided by 
placing the center in the insulating layer between two su- 
perconducting layers (first noticed by Bulaevskii^). The 
structure of the "core" is similar to the structure of the 
phase drop across a flux line in a two-dimensional Joseph- 
son junction^'^. This well-studied problem has a solution 
where the phase difference across the two layers drops 
by 2tt over a distance of the Josephson length Ajj^ For 
the 3D layered superconductor, this length is given by 
Aj = "fd, and we can think of a central region jd wide 
and d high as the core of an in-plane vortex. Beyond this 
core, the flux density and currents are quite similar to 
that for a continuous anisotropic superconductor^'*. The 
screening by z-axis currents is much weaker than that 
by in-plane currents, and the flux line is stretched into 
an ellipsoidal shape with a large width ~ Ac along the 
layers. Even though only the "core" resembles the vor- 
tex in a 2D Josephson junction, it has become common 
usage in the literature to label the entire flux line with 
this orientation a Josephson vortex. 

Consider now a flux line directed along the a;- axis. The 
general structure of this Josephson vortex was first de- 
scribed by Bulaevskii^: The center of the vortex lies 
between two layers, so that there is no core with sup- 
pressed amplitude of the order parameter, while at large 
distances from the center the structure is similar to a 
conventional flux line. The phase around the vortex is 
not given trivially by symmetry, but is a solution to the 



non-linear equations ^. The most convenient path to a 
quantitative solution is to separate the problem into two 
different scales: At large scales we can ignore the non- 
linearity and there is an analytical solution. At small 
scales the numerical solution is simplified by ignoring 
the screening contribution of the vector potential. Fortu- 
nately, for Xah/d ^ 1 there is a large region of intermedi- 
ate scales where both approximations work well, allowing 
us to match the small-scale and long-scale solutions. 

We consider a vortex centered between layers and 1, 
and J/ = 0, which is defined by the limiting values, 
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This corresponds to the following conditions for the in- 
terlayer phase difference 



^n,n+i = 0, for y -^ ±oo and n 7^ 0, 

0, for y — >■ +00, 
— 27r, for y — > —00. 
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To obtain the distributions of current and field, we first 
derive a useful exact equation for the magnetic field. The 
current components ^ and Q can be represented as 
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where 5"'"+^ is the average magnetic field between 
the layers n and n + 1 and V„ is a difference opera- 
tor, V„A„ = An+i ~ An- Collecting the combination 
(47r/c) {-Xl%Jn.n+i + {>ib/d)\7nJy,n), we obtain 
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with V„A„ 



^ ^ - An+i + An-i - 2 An. The difference of 
ipn+i,n and sin 93,1+1 _„ decays outside the nonlinear core 
and satisfies the relation 
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Therefore, in the continuum limit the right-hand side of 
([T3| converts into ^o^{y)^{z) and (fT3|) transforms into 
the usual equation for the vortex magnetic field^ 
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which gives 
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The current densities outside the core region are also 
given by standard formulas for anisotropic superconduc- 
tors. 
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These results should be valid as long as the linear ap- 
proximation for the sine of the phase difference is good. 
To find the range of applicability for this approximation 
we compare the last equation to ^, which gives near the 
vortex center. 
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(10) indicating that the linear theory breaks down at 



{y/ld) + n2 



1. This condition therefore sets the 



boundary of the nonlinear core. 

The above analysis shows that the Josephson vortex is 
characterized by two sets of length scales. A region where 
the interlayer phase difference is large defines the nonlin- 
ear core of the vortex. In the z direction this region is 
essentially localized within the central junction and in the 
y direction it spreads over the Josephson length "fd. At 
scales |z| , Ij/I /7 » d the vortex structure is described by 
the anisotropic London theory. In addition, in a wide re- 
gion one can neglect screening effects where the currents 
around the vortex decay as l/r, (though the current pat- 
tern is strongly stretched along the layers). Screening 
of the currents and magnetic field becomes important at 
the length scales l^l ~ Aab and \y\ w Ac, which are much 
larger than the corresponding boundaries of the nonlin- 
ear core. 

Due to this vortex structure, a quantitative analysis 
may be obtained with more ease by introducing an in- 
termediate scale i?int, with d < Rmt < Aab, such th at 
at a distance from the vortex center y^z'^ + (j//7)2 = 
^int both non-linearity and screening may be ignored. 
One then considers separately the small-distance region 
\/z-^ + (2//7)2 < i?int (containing the nonlinear core) and 
the large-distance region \/z^+{yfyP > Rint (where 
screening will become important). At small distances one 
can neglect screening. In the London gauge, V • A = 0, 
this means that the vector potential A can be dropped 
and the vortex is described in terms of in-plane phases 
4>n{y) only, which satisfy the following equation (from 
©) 

{ldf—r^'^s\n{(j)n+i " (/)„)- sin ((/)„ - (/)„_i) = (20) 
dy^ 

and the boundary conditions ([9]). These conditions are 
satisfied by our knowledge that outside the nonlinear 
core, {n — 1/2)2 _|_ (^yj^^y ^ -^^ ^j^g phase has to ap- 
proach the scaled version of the usual form relating to 



the angle around a vortex, 
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Multiplying (l20l) by dcjin/dy, summing over n, and per- 
forming an indefinite integral over y, we derive the fol- 
lowing exact relation for all y, 
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(22) 
which is analogous to the first integral of a second-order 
differential equation with one variable. For the case of 
an isolated Josephson vortex the constant is zero. In 
contrast to the single- variable case, this relation does not 
help us to find the exact solution of the coupled non- 
linear equations (|20|) . and one has either to use some 
approximate solution or to solve it numerically. Relation 
([^ can, however, be used to test the accuracy of the 
approximate and numerical solutions. 

A simple approximate solution has been proposed by 
Clem and Coffcy^^ (the CC solution) who used the fol- 
lowing ansatz for the magnetic field 
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and found that the best approximation for the core struc- 
ture is achieved by selecting the cut off j^cc = id/ 2. This 
field distribution allows one to obtain the distribution of 
the phase difference 
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In particular. 



this corresponds to (j^iiy) 



tan"^(7d/22/). 



The accurate numerical structure for the core was ob- 
tained in Ref. l56l . Figure[2]presents a visualization of this 
numerical solution, and we compare the phase difference 
in the central junction to that from the CC solution in 
figure [3] The numerical solution is characterized by the 
following properties. The maximum in-plane phase gra- 
dient is given by 
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2) and the 



(the CC solution gives ^d {A(f)i/dy)y=o 
maximum Josephson current fiows at a distance j/max — 
0.8470? from the vortex center (the CC solution gives 
2/max = ycc = 0.57(i). The maximum magnetic field in 
the vortex core is given by 
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FIG. 2. Visualization of the numerically computed structure 
of an isolated Josephson vortex. The arrows represent the 
current distribution (half interlayer distance corresponds to 
maximum Josephson current). The greylevel codes for the 
cosine of the interlayer phase difference. The scale in the 
y-direction is in units of the Josephson length Aj = 'yd. 
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FIG. 3. Sine of the phase difference between the central layers 
of the Josephson vortex. For comparison, the approximate 
solution of Clem and Coffey^^ is also shown. 



The asymptotic limits for the phase difference in the cen- 
tral junction are, 
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Outside the core one can calculate the correction 
(5(/)„ (y) to the continuum-limit phase asymptotics (|2ip by 
treating the discreteness and nonlinearity of the Joseph- 
son current perturbatively (see appendix IB)) . This gives 
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where R = \/{n — 1/2)^ -|- (y/jd)^, and the constant 
Cs^ ~ 4.362 is found from comparison with the numerical 
solution. 

We can find the energy per unit length of the Josephson 
vortex by inserting this solution into ([T]). The simplest 
method^^ is again to split the energy into two contribu- 
tions: one from the region at large distances where the 



linear approximation is valid, and one from small dis- 
tances where we need the numerical solution, but can 
ignore the contributions of A to the current (i.e. ignore 
screening). The first is found analytically, while the sec- 
ond needs a numerical integration. The final result is (see 
also^) , 



ejv = - [ln(Aab/d) + 1.55] , 

7 



(29) 



with £o — ^o/{4:'n'Xa,h)'^. This energy determines the 
lower critical field Hd^x above which Josephson vortices 
are generated. 



Hci,x = 47rejv/*o = 
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To summarize, the solution for a Josephson vortex 
presented here is very similar to the usual flux lines in 
isotropic superconductors, but stretched by the factor 7 
in the j/-direction. The reason for this similarity is that 
the linear approximation to the Josephson relation works 
well away from the vortex center. The important feature, 
however, is that at the center of the vortex there is no 
normal core, but rather a phase drop of nearly 2tt across 
the central junction over a distance of jd. 



A. Line-tension energy of Josephson vortex 

In this section we consider the line-tension energy of 
a distorted Josephson vortex, an important parameter 
which determines thermal wandering of the vortex line 
and its response to pinning centers. We consider a kink- 
free vortex located in between the layers and 1 and 
defined by the planar displacement field u{x). As the 
energy of the straight vortex does not depend on its ori- 
entation inside the layer's plane, for very smooth distor- 
tions with the wavelength larger that Ac the line-tension 
energy is simply determined by the line energy (j29p . 

6F = f dx^ (p^ for \du/dx\ < \u/Xc\ 

This simple result, however, has a limited interest, be- 
cause most properties of the vortex are determined by 
deformations with smaller wavelength, |du/(ix|/|u| ^ 
\kx\ ^ 1/Ac- In this range, the line-tension energy ac- 
quires nonlocality, a typical feature of vortex lines. An 
accurate calculation of the line tension for this regime 
presented in Appendix [K\ leads to the following result 
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with ej = Eo/^d and Ct ~ 2.86. The important feature 
is the logarithmic dependence of the efficient line tension 
on the deformation wave vector, which is a consequence 
of nonlocality. 



IV. DILUTE LATTICE, B^ < $o/2-!vyd^ 

When the Josephson vortices are well separated, the 
linear and continuous approximation can be applied to 
the energy functional ([Ij everywhere except in the core 
regions, which reduces it to the anisotropic London model 
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(32) 



This means that the lattice solution is just a linear ad- 
dition of single flux-line solutions and the lattice energy 
is determined by this London model. To understand the 
nature of the ground state it is useful to apply the rescal- 
ing trick^i^. 



(y,7z) and A = (Aj^,Aj7) 



(33) 



which in the case of zero z-component of the magnetic 
fleld precisely reduces the system to the isotropic state^i. 
Therefore the ground state configuration in scaled coor- 
dinates is given by a regular triangular lattice. In real 
coordinates this state corresponds to the triangular lat- 
tice strongly stretched along the direction of the layers. 
Within the anisotropic London model the lattice is de- 
generate with respect to rotation in scaled coordinates. 
In real coordinates this corresponds to an "elliptic rota- 
tion" illustrated in figure |4l In particular, there are two 
aligned configurations, in which Josephson vortices form 
vertical stacks along the z axis (see figure [5]). For these 
configurations the vertical distance between the Joseph- 
son vortices in the stacks, a^, and the separation between 
the stacks, Oj,, are given by 



a. = V/?$o/(7Sx), fly = V7*o/(/3S,) (34) 

where the constant /3 is equal to 2-\/3 and 2/-\/3 for the 
upper and lower configuration in figure [5] respectively. 

The interaction energy of the Josephson vortex lattice 
can be reduced to interaction energy of Abrikosov vortex 
lattice using scaling trick. This energy must be added to 
the self energy of each Josephson vortex ([25]) which, in 



the intermediate field regime Hc\,x ^ Bx 
the result 

^BJ BxEn ^^ A-23$o 
■'■^^^ Stt $o27 "^XjdWx 



<-B, 



,rf2, gives 



(35) 



The "elliptic rotation" degeneracy is eliminated by the 
layered structure of superconductor. There are several 
different mechanisms of this elimination. First, due to 
the strong intrinsic pinning, the vortex centers must be 
located in between the layers. This limits the possible 
lattice orientations. A second, less trivial, mechanism 
is from the corrections due to the discrete lattice struc- 
ture to the vortex interactions. The degeneracy is also 




(a) 



FIG. 4. Ground-state lattice configuration for an in-plane 
field and its rotational degeneracy within the anisotropic Lon- 
don model in (a) scaled coordinates and (b) real coordinates. 
The ellipse aspect ratio corresponds to the anisotropy factor 
~ 3, much smaller than the anisotropy of, e.g., BSCCO. 



Scaled coordinates 



Real coordinates 




FIG. 5. The two alternative lattice configurations that are 
aligned with the layers, in scaled and real coordinates. 



eliminated by thermal fluctuations, because the Joseph- 
son vortices mostly fluctuate along the layers directions 
and this selects preferential lattice orientations. All these 
mechanisms will be considered in details below. 



A. Selection of ground-state configurations by 
layered structure 

As the centers of the Josephson vortices must be lo- 
cated between the layers, the layered structure plays a 
crucial role in the selection of the ground-state lattice 
configurations. The Josephson-vortex lattice is commen- 
surate with the layered structure only at a discrete set of 
magnetic fields. Due to the "elliptic rotation" degener- 
acy of the lattice within the London approximation, the 
family of commensurate lattices includes lattices aligned 
with the layers (see Fig. [5]), as well as misaligned ones. 
To make a full classification of commensurate lattices we 
consider a general lattice as shown in Fig. [6l i^°i^^ . The 
lattice is characterized by three parameters: the in-plane 
period a, the distance between vortex rows in the z di- 
rection b = Nd, and the relative shift between the neigh- 
boring vortex rows in qa. The lattice shape is charac- 
terized by the two dimensionless parameters, q and the 
ratio r = b/a. The lattice parameters are related to the 
in-plane magnetic field, B^, as B^ = ^^/{ab). The two 
aligned structures in Fig. [5] correspond to (? = 1/2. As 
the replacement g -^- 1 — g corresponds to a mirror re- 




(m+n/2)ao 



FIG. 6. (a) General Josephson vortex lattice and its pa- 
rameters, (b) Orientation of layered structure with respect 
to ideal lattice (in scaled coordinates). The layered struc- 
ture fits the ideal lattice only if it is oriented along one of 
the crystallographic directions, which is characterized by two 
numbers {m,n), corresponding to expansion of the direction 
vector over the two basic lattice vectors, ei and 02. Several 
possible directions are shown with the corresponding indices 
{m,n). The layers, together with the lattice parameters, a, h, 
and q, are drawn here for the (3,1) orientation. 



flection with respect to the x-z plane, every structure 
with q ^ 1/2 is doubly degenerate. In addition to giving 
the general ground states, these lattices describe multiple 
metastable states with unique properties studied in Refs. 
rod and l6ll . which we will review below. 

We now classify the exactly commensurate lattices to 
give the set of commensurate fields. An equivalent ge- 
ometrical analysis has been done in Ref. i62l following a 
somewhat different line of reasoning, but with the same 
final result for the commensurate fields. The analysis 
of commensurability conditions can be done most conve- 
niently in scaled coordinates ([55]). In these coordinates 
the ground-state configuration corresponds to a regular 

triangular lattice with period oa — \/2j^o/\/3Bx. It 
is convenient to consider the orientation of the layered 
structure with respect to this lattice rather than the other 
way round. The layered structure fits this lattice only if 
it runs along one of the crystallographic directions, see 
Fig. Eh. This direction, {m,n), is defined by the lattice 
vector, em,n, which can be expanded over the two basic 
lattice vectors, e/„ „) = mei + ne2. For nonequivalent 
directions m and n must be relatively prime numbers 
(i.e., there is no integer other than one that divides into 
both m and n). Any such direction corresponds to a set 
of matching fields which we notate as i?(m,n)(A^)- We 
also notate the lattice parameters corresponding to such 
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an orientation as «(„„), b^rn,n}, and q(m^n)- Ininiediately, 
we obtain 



a(„i,n) = ei^^^n^ = ai~^\/vn?- + mn + \ 



(36) 



It is useful to write the unit vector perpendicular to 
the layers, z, in terms of e(m_„). This vector is labelled 
S(„_„) in Fig. Eb, and is given by, 



'(771, n) 



[7n.7i) 



X X 



-^{m.n) 



(37) 



Commensurability means that the projections of the two 
basic lattice vectors on S(„ „) should be an integer num- 
ber of layers, i.e., 

ei • S(„^„) = njd, 62 • S(„^„) = rhjd, (38) 

(in scaled coordinates the interlayer distance is 7^). Us- 
ing ([55)1 and ([57]) . we rewrite these conditions as 



V3. 



-OA 



2 y^m^ + TTm + n^ 

V3- 



-flA 



2 \/m'^ + mn + : 



— njd, 



= rhjd. 



(39) 
(40) 



These equations mean that rh/n — m/n. As m and n 
are by definition relatively prime numbers, the set of al- 
lowed m and n is simply given by to = Nm and h = Nn. 
Therefore we can represent the commensurability condi- 
tion as 

-— oa = N\/m^ + mn + ri?^d, (41) 

which gives the following set of commensurate fields, dis- 
tances between neighboring rows h = Nd and ratios 

^(m,n) 



V3 



$n 



Bira,n){N) ^ N^^d^[m^ + mu + u^Y 



^{m,n) 



— OA 



\/m'^ + mn + v? 



\/3/2 



771-^ -j- mn -{- n^ 



(42) 



(43) 



(44) 



Finding the parameter q[rn,n) for a general orientation 
is a more complicated problem. Defining the direction to 
the nearest-row site {mi,ni) (see figuretH]), we have 



H{m.n) ■ 



e{m.n)-e{rm,m) _ \mmi + (min+mni)/2 + nni\ 



■^{m^n) I 



m? -\- mn -t- n^ 



(45) 
Expressing the neighboring-row separation via (TO.i,r7i), 

|[e(m,n) xe(™i,„,)]| _ ^aA\min-mni\ 



^{va^n) 



^(m,n) 



^/w?' + mn + n^ 



and comparing it with Eq. (I43p . we can see that pair 
(7711,711) must satisfy the condition 



|77ii72 — mn\\ = 1. 



(46) 



It is well-known from the theory of numbers that for any 
relatively prime pair (771, n) one can find a complimen- 
tary pair (7771,77,1) satisfying this condition, and there is 
a general recipe to find complimentary pairs based on 
the Euclid algorithm (see, e.g., Ref. 63) . Moreover, as 
the combination 777177 — 77777,1 does not change with the 
substitution 7771 — >■ 7771 +m, ni — > 77i -I-ti, there is an infi- 
nite set of pairs which satisfy condition (1461) (physically, 
this corresponds to different lattice sites in the neighbor- 
ing row). Therefore, the problem to find qirn,n) can be 
formulated as follows: among all pairs (toi,77i) satisfy- 
ing condition (j46p one must find the pair which minimizes 
|777777i -|- (777177 -|- 7)777i)/2 -|- 7777i | and usc this pair in Eq. 
(|45|) . (Practically, we need not search to very high-order 
directions.) In the case 77 = 1 and arbitrary 777 the choice 
of (7771, 77i) is obvious, (7771,771) = (—1,0), and we obtain 



1(m,l) 



777 



1/2 



m^ 



(47) 



1? + 7?7 + 1 

We should stress that these results essentially rely 
on the linear London approximation, which implies 
a very strong inequality ca 3> 7^, or equivalently, 
N\/m?+rnn + n? 3> 1. The number of vortex- free lay- 
ers per unit cell is given by A^ — 1. The case A^ = 1 rep- 
resents a special situation when all the layers are filled 
with vortices and are equivalent. It is interesting to note 
that even for a dilute lattice one can have Josephson vor- 
tices in every layer (A^ = 1) in the case of high-order 
commensurability (777,77 ^ 1). In an ideal situation, the 
lattice transfers with changing magnetic field between 
different commensurate configurations via a series of first- 
order phase transitions. The number of competing states 
rapidly increases as the field decreases. 

A full analysis of the structural evolution requires con- 
sideration of the energy. In the London limit a very useful 
expression for the energy of the general lattice in Fig. [S^ 
has been derived in Ref. |60- We outline this derivation 
and present the final result in a somewhat different form. 
For the lattice in Fig. |B^ the interaction energy in the 
London limit is given by. 
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Jsing the formula 
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we can sum over k and integrate over y leading to 
i„t _ B^ 52 [TrAab sinh(6/Aab) 



/■111 
/jl 



8^ 2nXl^ 
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b cosh(6/Aab) - 1 
sinh [27r(7f,(/)] 



^ gb{l) cosh [2'Kgi,{l)] — cos {2Trql) 



dz—— 
9b{z) 



with gb{z) — y (5/27rAab) +r'^z'^ and r = b-f/a. This 
expression significantly simphfies in the intermediate re- 
gion b <C 27rAab, where we can use the expansion 



TrAab sinh (&/Aab) 
b cosh(6/Aab)- 



2^ 
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and drop b^ / (27rAab) hi gb{z) meaning that gb{z) — > rz. 
This allows us to represent the interaction energy in this 
regime as^ 
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2 \27rAabAci?2:, 

with 7e = 0.5772 being the Euler constant and 



(48) 



^ , . TTT -r-^ cos(27rg/)-exp(-27rrZ) 1 , ,„ ^ 
gL(r,g) = - + )_^ ^[cosh(27rr0-cos(27rg0] " 2 ^"^'""^- 

(49) 
The dimensionless function Gl {r, q) depends only on 
the lattice shape. Its absolute minimum corresponding 
to the triangular lattice is given by Gl(-\/3/2, 1/2) ~ 
—0.4022. A peculiar property of the function GL{r,q), 
following from the rotational degeneracy, is that it 
also has this value for the whole set of pairs (r, q) = 
(j'{m.n)i<l{m,n)) Corresponding to the different lattice ori- 
entations. In particular, for (m, n) = [m, 1) we have 



G 



73/2 



ni+1/2 \ _ 



2^ ^1, 2a. ^1 / - Gl(V3/2,1/2). This func- 
tion also has very peculiar behavior at small r which is 
important for the statistics of metastable states*'^: at 
r — > it acquires peaks at all rational values oi q — k/l. 
Large-order peaks with denominator / develop as r drops 
below l/(27r/). 

For layered superconductors we have 



b = Nd, a = 
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with B^ii2 = $0/(7'^^) and, adding the energy of isolated 
Josephson vortices, we can write the total energy of the 
lattice as 

f,,{N,q,h) = ^ 
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1.432 + GL(r,(j) 



(50) 



with h = 2ttBx/Bj^2 and r ~ N'^h/2n. For given h 
the ground state configuration is determined by the min- 
imum of Gi,{N'^h/2Tr,q) with respect to discrete N and 
continuous q. As follows from Eq. (|42|) . perfect fits where 
Gl reaches its absolute minimum occur at the set of re- 



duced fields h 
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(51) 



At these fields this energy reproduces the result ((35)) . The 
field dependence of Gl for the ground state is shown in 
Fig. m The continuous London model does not not ac- 
curately describe layered superconductors at high fields. 
To obtain lattice structures in this region one has to con- 
sider the more general Lawrence-Doniach model. The 
transition between the aligned lattices have been studied 
within this model by Ichioka^l. However, our analysis 
in the next section shows that at many fields the true 
ground state is not given by an aligned lattice. 



Evolution of ground-state configurations within 
Lawrence-Doniach model 



The accurate analysis of the lattice configurations 
within the Lawrence-Doniach model which we report in 
this section was only published in short ProceedingSS.. In- 
dependently, such numerical analysis was done by Nono- 
mura and Hu^^ with fully consistent results. 

At high in-plane magnetic fields the spatial variations 
of the field are very small and in the first approxima- 
tion can be neglected. In this limit the only relevant 
degrees of freedom are the superconducting phases and 
the relevant part of the LLD energy ([1]) per unit volume, 
/^ = Fi^i,u/{LxLyLz) - Bl/iSir), can be written as 



U[Mr)]^^^Y.Jdy 
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(52) 



To simplify the analysis we will introduce the reduced 
in-plane length y = y/^d and the reduced magnetic field 
h = 2'K^d'^Bx/^Q leading to 



UVt^nir)] = 



1— cos ( 
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2 (Vj,(/)„)^ 



(53) 



with £j = Eg/^d. Varying this energy with respect to the 
phases (t>n{y) we obtain an equation for the equilibrium 
phase distribution [equivalent to ([2]) when we ignore the 
spatial dependence in B^] 



^h 



-sin (0„+i -(t)n + hy) - sin (0„ - (/)„_i + hy) = 0. 

(54) 

We again consider a general lattice shown in Fig. [6]i 
with in-plane period a, N layers between neighboring 
rows, and relative shift qa between relative rows with a 
and iV being related to the reduced field as /i = 2'iT^d/Na. 
It is sufficient to find the solution for the phase in one 
unit cell, 0<2/<a,l<n<A^, using appropriate quasi- 
periodicity conditions for the phase. The total lattice 
energy per unit volume can be represented as 



S,$o 



(4^)2AabA, 



■u{N,q,h), 



(55) 
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where the reduced energy u{N, q, h) per unit ceh is given 

by 



uiN,q,h) 



-1 — cos ( 



1 ^ r 



n=l 



dy 
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hy)] . 



(56) 



Using a relaxation method to solve (l54l) numerically 
within one unit cell, we can find the energy u for any 
given values oi N, q and h. To match the London repre- 
sentation (|151) . we write u{N, q, h) in the form 



uiN,q,h) 



2 h 



IA323 + G{N,q,h) 



(57) 



where the function G{N, q, h) defined by this equation 
approaches the London limit Gi,{r = N'^h/2Tr,q) for h — > 
0. 

We consider first influence of the layered structure 
at small fields. As shown in the appendix [B] in the 
lowest order with respect to h the layered structure 
gives orientation-independent correction to energy, G « 
{h/32)ln{Ch/h). In higher (quadratic) order the layered 
structure generates orientation-dependent correction to 
the lattice energy leading to a break down of the "elliptic- 
rotation" degeneracy of the lattice. To study this ef- 
fect quantitatively we plotted in Fig. [7] the computed 
field dependencies of G{N, q, h) for several lattice ori- 
entations at the corresponding reduced commensurate 
fields /i(m,„)(A^) given by (I5T|). At smaU h, h < 0.05, 
neglecting a very weak dependence on orientation, one 
can accurately fit the correction from the layeredness as 

G(/i) - Gl w ^ In ^ - 0.075/1^. (58) 

One can see that among the two aligned structures shown 
in Fig. [5] the layers favor the lower structure with indices 
(1,1). However, for ft, < 0.1 the energy difference between 
the two structures is tiny and in real samples external 
factors may select the lattice orientation. On the other 
hand, one can expect that at sufficiently large fields the 
ground-state configuration will be selected by the layered 
structure even in real samples. 

Energy corrections due to the layered structure fa- 
vor lattice stretching along the layer direction and shift 
down the matching fields. This effect is strongest for the 
aligned lattice (1,0) and is illustrated in Fig. |S1 In this 
figure we show the field dependences of G{N, h, 0.5) for 
different N. When a smooth function is subtracted from 
these dependences, local minima are realized at fields 
h{ifi){N) which are smaller then the London matching 

field /i(i^o)(^)- The shift ft(i^o)(^) - ^(i,o)(^) rapidly 
decreases with increasing magnetic field. We found that 
this shift is described by the following equation 



1(1,0 



~{N) 



huMN) 



For other lattice orientations the shift is smaller but still 
noticeable. To quantify the energy difference between 
the aligned lattices due to the layered structure, we fit 
their energies at the shifted matching field for /i < 0.1 to 
smooth curves and subtract these curves. This procedure 



gives G(i,i) - G(i 



^(1,0) 



-0.011^2. 



l + (0.63/iV2)ln(19//i(i,o)(A^)) 



We can now explore the evolution of the ground- 
state configuration by direct minimization of the en- 
ergy with respect to the lattice parameters N and q 
as defined in Fig. [5] To this end we have computed 
the reduced ground-state energy defined as G{h) = 
min ff^q[G{N,q, h)]. We checked that if we consider only 
aligned lattices, the results of IchiokaM are reproduced 
for the transition fields between lattices with different 
periods iV in the case of large anisotropy. For compar- 
ison, we also made a similar calculation for the London 
model and computed the field dependence of the function 
Gl(/i) = minAr,,[GL(r = N^ h / {2tt) , q)] where Gi^{r,q) 
is defined in Eqs. p8)) and (|49)) . In figure [9] we com- 
pare field evolutions of these ground-state reduced ener- 
gies and the corresponding c-axis period N. For clearer 
comparison we subtracted from G{h) its fitted correction 
from Gl(\/3/2, 1/2) at small h given in ([55]) . Values of 
the London commensurate fields ft(m.ri)(-^) are shown on 
the top axis with several low-order fields being marked 
by corresponding indices using the format {m,n)pf. As 
expected, Gi^{h) reaches its absolute minimum for every 
h(m,n) {N) ■ We can observe several interesting properties. 
As the lattice orientation with indices {m,n) = (1,0) is 
not favored by the layered structure, several low-A^ con- 
figurations, 3 < A^ < 6, expected a.t h — /i(i o)(-^): are 
skipped. However, as one can see from the inset in Fig. 
[TUl for N — 5 and 6 the ground-state energy is smaller 
than the energies of these states a.t h — /i(i o)(^) only 
by a tiny value. For ft < 0.2 the actual evolution of lat- 
tice structure starts to follow roughly the London route 
(except for skipped state (1,0)6 near ft = 0.16) but with 
small negative offset, i.e., we again see that the match- 
ing fields systematically shifted down in comparison with 
their London values. 

The field dependence of the energy function G{N, q, ft) 
in an extended field range is shown in figure [TU] for the 
ground state and competing states. Each curve corre- 
sponds to the minimum of G{N^ q, ft) with respect to q 
at fixed ft and N and is marked by its value of N. We also 
show the first six lattice configurations which are realized 
with decreasing field. The inset of the figure blows up the 
low-field region. One can see that many lattice configu- 
rations compete for the ground state at small fields and 
at several fields (e.g., at ft w, 0.19, 0.137, 0.105 .. .) one 
or more lattice configurations have energies very close 
to the ground-state energy. We also note that there 
are several extended field ranges where in the ground 
state all layers are homogeneously filled with vortices 
{N — 1) even in the region of the dilute vortex lattice, 
e.g., 0.115 < ft < 0.17, 0.21 < ft < 0.38. 

We see that an accurate consideration within both 
London and Lawrence-Doniach models shows that the 
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FIG. 7. -Left panel: Field dependence of the reduced-energy function G(N , h, q) for several lattice orientations [m, n) at the 
commensurate field /i(m,n)(-^)- To enlarge small differences, we plot in the right panel the difference between G and its fit 
obtained using all data for h < 0.05. 





-0.32 




-0.34 


o 






-0.36 




-0.38 




-0.40 




0.000 






o 


-0.002 


xf 


-0.004 






CD 


-0.006 






(- 






-0.008 


CD 











FIG. 8. Upper panel: The field dependences of the reduced- 
energy function G = G{N,h, 0.5) for the aligned lattice (1,0) 
and different A'^. Vertical bars mark locations of the London- 
model commensurate fields /i(i^o)(A/"). Lower panel shows dif- 
ference between G{N, h, 0.5) and smooth curve through the 
points (/i(i,o) (■^)i G{N, /i(i,o) (■^)i 0.5)) (dashed line in the up- 
per panel). One can see that the matching fields systemat- 
ically displaced to the lower values h^iQ^{N), as illustrated 
for N = 4. Inset in the upper panel shows lattice structure 
at the displaced matching field for N = 3 (solid symbols) in 
comparison with the regular-hexagon structure at the London 
matching field. 
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FIG. 9. Upper panel: The field dependence of the reduced 
energy functions for London model (Gl, upper curve) and full 
Lawrence-Doniach model (G — Gat, lower curve). For clearer 
comparison we subtracted from G(h) its fit at small h given 
by Eq. (I58|l . Values of commensurate fields /i(m,„)(A'^) are 
shown in the top axis and the corresponding indices for several 
of them are written in the format {m,n)M- As expected, 
Gl reaches its absolute minimum for every hi^^^„){N). The 
lower panel shows field dependence of A'^ for ground state for 
both models (stripes for the London model and circles for 
the Lawrence-Doniach model). The same graylevel codes the 
value of A'' in the upper panel and the London-model plot in 
the lower panel. 



ground state of the Josephson vortex lattice at low 
temperatures does not give any preference to the lat- 
tices aligned with the layers. Therefore for equilib- 
rium field dependencies one can not expect to observe 



any strong features at the matching fields of these lat- 
tices, S(i^o)(7V) and B(i_i)(iV) given by Eq. ^. Nev- 
ertheless, clear commensurability oscillations have been 
observed experimentally in underdoped YBCO in irre- 
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FIG. 10. Field dependence of the energy function G(Af, g, h) 
for the ground state and competing states. Each curve cor- 
responds to the minimum of G(A^, q, h) with respect to q at 
fixed h and A''. The curves are marked by the value of A^. 
Lattice configurations in scaled coordinates are shown at six 
marked fields. Inset illustrates competition between different 
configuration at smaller fields 




FIG. 11. Levitov's hierarchical plot of metastable states 
in q-r plane— (in this plot q is selected within the inter- 
val [0.5,1]). Each dotted curve is obtained from the lo- 
cal minima of function Gl {r, q) with respect to q at fixed 
r. New branches appear as a result of "quasi-bifurcations" . 
Each "quasi-bifurcation" is associated with a rational number. 
The branches turn at points {q(m,n),T{m,n)) corresponding to 
ground states (marked by squares and labelled by the indices 
mn in the plot). 



versible magnetization^^''^'^ and nonlinear resistivity''^. 
The period of these oscillations corresponds to the fields 
B(ifi){N) indicating that for some reason in this mate- 
rial the aligned lattice (1,0) occurs to be preferable. We 
note that, due to small differences between the energies of 
different configurations, in real materials aligned lattices 
can be selected by external factors, such as interaction 
with correlated disorder (twin boundaries or dislocations) 
or sample surface. We will also see that the aligned lat- 
tice with indices (1, 0) is favored by thermal fluctuations. 
Finally, we mention the work of Ikeda and IsotaniS who 
performed similar analysis of the ground state configura- 
tions for field applied along the layers within the lowest 
Landau level approximation. 



Properties of metastable states in London 
model 



Josephson vortices can slide easily along the layers but 
there is a huge barrier for the motion across the lay- 
ers. This property makes it hard to equilibrate the lat- 
tice. It also leads to the appearance of a very large num- 
ber of metastable states. The properties of these states 
have been considered in Refs. (GOi and 61. Systematically, 
metastable states at fixed c-axis period can be sampled 
by first slowly cooling down the superconductor at fixed 
magnetic field and then in a second step decreasing the 
magnetic field at a low temperature—. We assume that 
the prepared starting configuration is the aligned lat- 



tice. As the c-axis period, N, is locked by the layers, 
the lattice stretches along the layers with lowering the 
field, i.e, the ratio r = b/a decreases. During stretching, 
these fixed-A*" metastable states go through a sequence 
of nontrivial structural transformations. In the London 
regime, the aligned configuration becomes unstable at 
To « 1.51/(27r) « 0.24*^0. This instability is driven by 
the repulsion between neighboring vortices in the verti- 
cal stack. At low r < To the parameter q continuously 
decreases starting from 1/2 to lower values. We found 
that the layeredness stabilizes the aligned structures: the 
critical ratio decreases to 0.231 at A*" = 3 and to 0.224 at 
N = 2. It is important to note that the shear instability 
occurs in the ground state only for A'^ = 1 (we consider in 
detail this structural phase transition below). At higher 
values of A^ this instability only occurs when the state for 
this given A^ is metastable with respect to other values 
of A^. This instability is considered in detail below. 

The statistics of metastable states has been explored 
in detail in Ref. [6l| where the similarity to the phyllotaxis 
phenomenon in biological systems has been pointed out. 
For every r one can find all local minima, qi{r) of the 
energy function G{r, q) with respect to q and plot all 
these minima in the q-r plane (see Fig. Ilip . The ob- 
tained pattern is quite peculiar. At r > rg the only 
minimum is at qo{r) = 1/2. Below r = rg this trajectory 
symmetrically splits into two. With further decrease of 
r many more minima appear forming a complex hierar- 
chical structure. The pattern can be viewed as a series 
of "quasi-bifurcations" occurring near rational values of 
q. "Quasi-bifurcation" corresponds to the appearance of 
a new branch below a certain value of r in the vicin- 
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ity of the old branch. The branches turn at the points 
('Z(m,n)i''(m,n)) Corresponding to ground states. The evo- 
hition of the initial state is described by the two main 
trajectories symmetrically split from q = 1/2. The tra- 
jectory with q > 1/2 "quasi-bifurcates" a,t q — Fj/Fj+i 
where Fj are the Fibonacci numbers and approaches the 
"golden ratio" (a/5 - l)/2 « 0.618 as r ^- 0. It goes 
through ground states with the indices also described by 
the Fibonacci sequence, {m,n) = (Fj^i, Fj). Unfortu- 
nately, these exciting predictions have never been verified 
experimentally because there is no direct way to probe 
the structure of the Josephson vortex lattice. 



D. Elasticity of dilute Josephson vortex lattice 

Josephson vortices easily slide along the layers but mo- 
tion across the layers is strongly suppressed by intrinsic 
pinning from the layers. Due to the intrinsic pinning, z 
axis fiuctuations of the vortex lines occur via kink for- 
mation. In moderately anisotropic layered superconduc- 
tors, such as YBCO, in which the c-axis coherence length 
is larger than or comparable with the interlayer spac- 
ing d, the intrinsic pinning potential V{uz) can be de- 
scribed as a cosine function of z axis vortex displacements 
V{uz) = ~Vo cos{2ttu z{x)/d). However such description 
becomes inadequate in strongly layered materials, where 
the structure of a kinks is very similar to the structure 
of a pancake vortex. 

In strongly layered materials at low temperatures 
one can neglect kink formation and take into account 
only in-plane lattice deformations u(r) = Uy(r) (planar- 
fluctuations model). In this case, one can derive the fol- 
lowing nonlocal elastic energy in the fc-space 

1 f dPii 

J^e/ = 2 y (^ [Cll(k)fc^ + C44(k)fc2 ^ ^gg^2] |^(]^)|2 

with elastic moduli 



C66 



Cll(k) = 
C44(k) = 



(87r)2Ah' 



S,$o 



l + Xl^k^ + XU^l + kl) (8^)2AabAc 
l + Xl^k^z + XHkl + kl) 



B^% 



■In 



(^"^'^^'^^^ s^c^' + hk^/nf 



(59) 

(60) 
(61) 

(62) 



While the tilt [c44(k)] and compression [cii(k)] moduli 
are not sensitive to exact lattice structure, the formula for 
the shear modulus cgg is valid only for perfect matching 
between the Josephson vortex lattice and layered struc- 
ture, which is achieved at matching fields (H^ . For a 
general lattice shown in Fig. |B^ one can derive a more 
general expression for cge using representation (|48l) " (|49p 



for the lattice energy^ and relation between lattice de- 
formation and change of parameter q, 6q — rdu/dz, 



C66 



B,^o , . 

s2 .n 966{r,q) 



iSnrXh" 



(63) 



with 



geeir,q)=4r^^GLir,q) =-{47Tfr'^ 



E 

1=1 



dq^ 
cos{2'Kql) cosh(27rrZ) — sin^(27rg/) — 1 



(cosh(27rr^) — cos(27r(7^)) 



/ sinh(27rrZ). 



This formula reproduces the result ()60|) for the commen- 
surate configurations (r, q) = (?'(„,,„), (/(n.m))- It also de- 
scribes instability of aligned configuration {q = 1/2) at 
r « 0.2^. 

Softest mode in the planar model corresponds to shear- 
ing between neighboring planar arrays of Josephson vor- 
tices. The harmonic approximation breaks for this mode 
first. The simplest extension of the linear elastic energy 
which describes strong interplanar fluctuation is amounts 
to replacement of the continuous displacement field u(r) 
by displacement of the planar arrays Uj {x,y) = u{x,y, jb) 
and replacement the shear term in the energy by the non- 
linear interaction term 



C66 



"■•'^ I T.) 



du\ ' 



— > 



2 
,2 " ^66 

d r- 



(27r)' b 



E 



1 — cos 2tt 



ij+l-Uj 



Such extension has been used to study strong- 
fiuctuations region^. 



V. DENSE LATTICE, B^ > '!>o/2TT'yd^ 

With increasing magnetic field distance between 
Josephson vortices decreases and at field B ^ B„ — 
^o/2TTjd'^ = BjcI^/2tt this distance becomes of the order 
of the vortex-core size. In contrast to the Abrikosov vor- 
tex lattice, for which overlap of the vortex cores marks 
disappearance of superconductivity, for the Josephson 
vortex lattice this field just marks a crossover to a new 
regime, the dense Josephson vortex lattice. Existence of 
this regime was pointed out by Bulaevskii and Clem'^^. 
In the dense Josephson vortex lattice the gauge invariant 
phase difference is a smoothly increasing function of dis- 
tance and the Josephson coupling energy can be treated 
as a small perturbation. This allows for the following 
quantitative description. 



A. Very high fields: Quantitative description using 
expansion in Josephson coupling 

At high fields Bj; > B^ , vortices homogeneously fill all 
of the layers. This means that all layers are equivalent 
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and the in-plane lattice period is a = 27r//i (see figure [T2l) . 
When the strong inequality B^^ » -Bcr {h >> 1) is satisfied 
Eq. (|54l) for the phases can be solved using an expansion 
with respect to the Josephson currents. In the zeroth 
order we can construct a regular lattice with an arbitrary 
translation from layer to layer by using the form, 

This corresponds to the gauge-invariant phase difference 

<pi°L+i = Kn + hy, 

i.e., the planar lattices in the neighboring layers are 
shifted by the fraction q — k/2'k of the in-plane lattice 
spacing a. In the first order we obtain 

^I'/'n ^ + sin {nn + hy) - sin (^(n - I) + hy) = Q 

which gives 

</"« '' (y) = T^ [sin {i^n + hy) - sin (^(n - 1) + hy)] . 
h-' 

Substituting this solution into ([55]) . we obtain the energy 
per unit volume up to second order with respect to the 
Josephson coupling, 



Ui^.h) 



-f<P 



T / 1 — cos K 

^ ' l- 



2/^2 



(64) 



We can immediately see that the minimum energy 
fmin{h) = (£j/7(i^)(l — l//i^) is achicvcd at k = tt, cor- 
responding to the triangular lattice shown in figure 1121 
The phase distribution in the ground state is given by 

0n(2/) « Trfc^ + ^Hl: sin {hy) (65) 

From this solution we can recover the distributions of the 
in-plane and Josephson current 



JvAv)'' 

jz,n{y)'' 



1)" /2TTdB^y 

7JJ cos 



h 
-l)"jjsin 



$n 



4(-l)" . f2TrdB^y 
■ sin 



2T:dBxy 



h^ ^"'\ $0 
and a weak modulation of the in-plane field 

' 2TrdBxy 



$n 



B,{y)^B, 



1)"$§ 



S,(2^dAe)2 



$n 



A schematic distribution of the currents is shown in fig- 
ure [HI 



B. Dense lattice close to the crossover region. 
Structural phase transition 

When the magnetic field approaches the crossover field 
^o/{2TT^d'^), the perturbative approach of the previous 
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FIG. 12. Schematic distribution of currents in the dense 
Josephson vortex lattice. The circles marlc the centers of the 
Josephson vortices. 



section becomes insufficient and one has to obtain a full 
solution of the nonlinear equation (|54p . The general so- 
lution for the lattice with the arbitrary phase shift k can 
be written as 

, / X n(n — 1) / Kn\ 

(^n{y)^n "■ ^ > +g[y+—y (66) 

where g{y) is a periodic function, g(jj + 2TT/h) — g{y), 
which obeys the following equation 



d^g 



sm 



g (y) + hyj 

{g [y) -g(y-^+hy-n)^o. (67) 




The reduced energy / = f^^d^ /ej can also be written in 
terms of the function g{y), 



(68) 



Equation (|67p does not have an analytical solution and 
has to be solved numerically. Lattice configurations of 
the dense lattice also has been investigated using the code 
developed for the lattice with general period N. Both ap- 
proaches give identical results. Numerical investigation 
shows that the triangular lattice with k = n gives the 
ground state for h > 1.332. At /i sa 1.332 the system has 
a second-order phase transition to a lattice of lower sym- 
metry, see lattice structures for h = 1.35 (a) and h — 1.2 
(b) in figure [TUl The field dependence of k and corre- 
sponding lattice shift q are shown in figure 1131 Ikeda 
and Isotani^^ found that within the lowest Landau level 
approximation this structural phase transition occurs at 
somewhat higher value, /i w 1.4. 

To study validity range of the high-/i approximation of 
the previous section we plot in figure [T3] the computed 
field dependence of the reduced energy together with its 
high-field asymptotics, derived in the previous section. 
As one can see, the perturbative approach gives a good 
approximation for the energy down to h ^ 2. 
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FIG. 13. Field dependence of the phase shift k and corre- 
sponding lattice shift q for the dense Josephson vortex lattice. 
At h ~ 1.332 the lattice experiences a continuous structural 
phase transition. 
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FIG. 14. Field dependence of the reduced energy for the 
dense Josephson vortex lattice. The dashed line shows the 
high-field asymptotics. The arrow marks the position of the 
structural phase transition at /i ~ 1.332. 



C. Elasticity of dense lattice 

In this section we consider the deformation energy 
of the dense Josephson vortex lattice in the limit h ~ 
2ttj(PBx/^q ^ 1. In particular, this energy serves as a 
starting point for analysis of fluctuations. We will fol- 
low the approach used by Korshunov and Larkin®^. The 
starting point of the analysis is again reduced LLD en- 
ergy in the phase approximation (|53p which we rewrite 
now for the general case of the phase depending on both 
reduced coordinates f = (x, y) = r/jd, 



n 



d^V 



T;\—7^] ~COs{(l}n+l~4'n + hy) 



(69) 



The ground-state phase distribution is given by Eq. (1651) . 
Now we consider small deformations of the lattice and 
split the total phase into the smooth (i;„) and rapidly 
oscillating in the y direction {4>n) parts 



(t>n{r) 



n{n + 1) 



-w„(r) + (^„(r), 



(70) 



where we assume dvn/dy <^ w„ and 0„ <Si 1. As 
the smooth part of the gauge-invariant phase difference 
is given by h{y + (u„+i — f„) /h) + im, the quantity 
M„ = — {vn+i — Vn) /h represents a local lattice dis- 
placement. Substituting representation (j70p in the en- 



ergy (|69)) . expanding with respect to 0„, and dropping 
rapidly-oscillating terms, we obtain 



Fj^KiEi 



E 



dr 




+ ((t>n+l-<i>n) Sm{Vn+l~Vn + hy + TTn) 



(71) 



As (J3n rapidly oscillates only in y direction, we kept only 
its y derivative. Minimization of this energy with respect 
to 0„ gives 

+ hy) + sin (vn-Vn-i + hy) 



K-1) 



„sin>^ 



n+l' 



/l2 



Substituting this solution into Eq. ([7T|) and averaging 
with respect to y, we finally obtain the coarse-grained 
energy of the deformed dense Josephson vortex lattice^, 
which we write in real units 



Ft 



4i: 



dr 



dv„ 
dr 



cos {v„ 



f^n-l-l" 



'2w„)-H 



(Aj/i)2 



(72) 

This energy describes the phase fluctuations in large in- 
plane magnetic field. The first term is just usual in- 
plane phase stiffness energy. In the elasticity-theory lan- 
guage this term represents the compression (dvn/dy) and 
tilt [dvn/dx) contributions. The second term represents 
the shearing interactions between the Josephson vortex 
arrays in neighboring junctions. It originates from the 
Josephson coupling energy and can be viewed as the ef- 
fective Josephson coupling renormalized by the in-plane 
magnetic field. Roughly, we can state that with in- 
creasing magnetic field the effective Josephson energy de- 
creases as l/K^ and the effective Josephson length, Aj^i, 
increases linearly with h, 



Kjh = Ajh = 



2nj^d^B^ 

$0 



(73) 



In the case of slowly-changing from layer to layer defor- 
mation, we can expand cosine in Eq. (|72p and obtain the 
harmonic elastic energy of the dense Josephson vortex 
lattice in terms of smooth phase deformations 



F, 



(f>—el ' 
Eq 

2d 



2 ^ 



^> rfr 



dVn 

"d7 



[2^f 



''1^ dk^ 
-ttM 2tt 



kt 



2(A.7/l)2 

2 (1 — cosfczd)^ 



{Ajhf 



,fc|2 



(74) 



(75) 



Using relation between the phase perturbation and lat- 
tice displacements 

k '*■"*' 

V — 

Aj (exp(zA:zd) - 1)' 

we can rewrite the elastic energy in a more traditional 
way, via lattice deformations 

1 /• d^kii r''^ dk^ 
2] (27r)2./ ^,^^ 



F^ 



ei ^ 77 



TT/d 



ciiik,)kf + ceeki lu"] 



,fe|2 



(76) 



16 



1. Rectangular lattice 




FIG. 15. Size-magnetic field phase diagram of tlie confined 
Joseplison-junction stacli. Daslied line separates the large- 
size and small-size regimes. Black lines correspond to inte- 
ger flux quanta per junction. Shaded areas mark regions of 
rectangular-lattice ground state. Representative lattice con- 
figurations in two points are illustrated by plots of oscillating 
Josephson currents in two neighboring layers. Small ellipses 
mark the centers of the Josephson vortices. 



with elastic constants 
Bl 



cii{k^) 



4^ klK^. 



C66 



^l 



i2T:^(P-1^\l^ 



where we used notation kz = 2sin(fczd/2) /d. Note that 
in our case the nonlocal tilt modulus Ci4{kz) is identical 
to compression modulus cii(fcz) and they coincide with 
elastic moduli within the anisotropic London model (pT|) 
and (j62p in the limit fc^Aab 3> l,fca;A^. This elastic en- 
ergies ([711) . (1751 . and (1751) can be used to study weak 
fluctuations and weak pinning of the dense Josephson 
vortex lattice. The shear modulus is field independent in 
the dense-lattice regime. One can check that it matches 
the dilute-lattice result (|60l) at the crossover field. 



D. Lattice configurations and magnetic oscillations 
in finite-size samples 

In this section we consider dense-lattice configura- 
tions in finite-size samples. This study is actually mo- 
tivated by experimental observations of magnetic oscil- 
lations in small-size BSCCO mesas with lateral sizes 2- 
20 /im^"— . Such small-size mesas behave as stacks of 
intrinsic Josephson junctions with strong inductive cou- 
pling between the neighboring junctions. The detailed 
analytical theory describing the magnetic field depen- 
dences of lattice configurations and critical current has 
been developed in Refs. [TO. Lattice structures also have 
been extensively explored numerically^^^i^i^ and both 
approaches give identical results. In a small-size sam- 
ple the lattice structure is determined by two compet- 
ing interactions: the interaction with boundaries that 
favors an aligned rectangular configuration and the bulk 
shearing interaction between neighboring layers which fa- 
vors a triangular configuration. Depending on the mesa 
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FIG. 16. Illustration of the oscillating magnetic field depen- 
dence of the critical current for L — 4Aj. Crossover between 
$0/2 and $0 periodicity is seen at hL = B^/Bl ~ 1. Shaded 
areas show the regions of stable rectangular lattice. 



width L and magnetic field, two very different regimes 
realize. In the large-size regime the vortex lattice is tri- 
angular and it is only deformed near the edges. In the 
small-size regime the lattice structure experiences a peri- 
odic series of phase transitions between rectangular and 
triangular configurations. The triangular configurations 
in this regime are realized only in narrow regions near 
magnetic-field values corresponding to an integer number 
of fiux quanta per junction where the interaction with 
edges vanishes. The typical width of the mesa which 
separates these two regimes is given by the length Ajh, 
Eq. (|73|) . which is proportional to the applied magnetic 
field. Hence, the crossover from one regime to another 
is driven by the magnetic field and the corresponding 
crossover field scale is Bl = Ba-L/Aj = L$o/(27r7^d^), 
for B^ > Bl the small-size regime realizes. The size- 
field phase diagram is shown in Fig. 1151 The regimes are 
characterized by distinctly different oscillating behavior 
of the critical current as function of the magnetic field. In 
the small-size regime, the critical current oscillates with 
the period of one fiux quantum per junction, similar to 
a single junction. In the large-size regime, due to the 
triangular lattice ground state, the oscillation period is 
half fiux quantum per junction. 

The quantitative study of the described behavior is 
based on the reduced energy ([M|) . which have to be 
rewritten for the finite-size case, < y < L = Ly, 
and also assuming that the system is uniform along the 

field direction, i.e., Jdr — )• L^ J^ dy. This energy has 
to be supplemented with the boundary conditions at the 
edges, d(j)n/dy = for y = Q,L. The important param- 
eter in the case of finite-size sample is the total mag- 
netic flux through one junction, $ — B^dL, which is 
connected with the reduced magnetic field by relation 
hL — 27r$/$o- In the dense-lattice limit, we again use 
presentation of Eq. ((70)) containing the smooth phase Vn, 
and rapidly oscillating component 0„ . It is natural to as- 
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sume that the interactions with the boundaries preserve 
the alternating nature of the vortex lattice. In this case 
symmetry allows us to take the smooth phase in the form 



Vniy) = an + {-iy'v{y) 



(77) 



where a describes the translational displacement of the 
lattice and v describes lattice deformations with respect 
to the triangular lattice. In particular, one can show that 
the maximum value of v{y), Wmax = 7i'/4, describes the 
rectangular lattice, i.e., identical (/>„ in all layers up to 
27r phase shift. The corresponding rapid phase becomes 
(j)n{y) ~ (— 1)"2 cos(2u) sin (a +hy) /h? . Averaging with 
respect to the rapid oscillations for such Vn{y), gives the 
reduced energy per layer and per unit length along x, 
U^F^Aj/iNL^Eo), 



— — [sin (2110) cosa — sin(2ui)cos (hL + a)] 
h 



dy 



1 + cos(4w) 



/l2 



(78) 



in which the bulk part directly follows from Eq. ([7^ for 
general Vn(y)- Variating this energy with respect to v{y), 
we obtain that it obeys the static sine-Gordon equation 










(79) 



with the boundary conditions 

dv 2 

—^(0) = — — cos(2wo)cosa, 

dy h 

-rz{L) = — — cos(2wi)cos ihL + a) 
dy h 



(80) 



Substituting solution of these equations into the energy 
functional (175)) gives the energy as a function of the lat- 
tice shift a, f^{a). Minimum of the energy with respect 
to a gives the ground state for given h and L. Higher- 
energy states at other values of a typically carry a finite 
current. The total Josephson current flowing through the 
stack is proportional to df^/da. Taking derivative of the 
functional (j72p with respect to a, assuming that at ev- 
ery a it is minimized with respect to v{u), we obtain the 
total current in units of jjAjL^ 

J{a)^— [sin(2vo)sina— sin(2ui) sin(/iZ-|-a)] . (81) 

An important consequence of this equation is that 
nonzero current exists only if the surface deformations 
vq and vl are finite. 

The general solution of equations ([7^ and ([M]) can be 
written in terms of the elliptic integrals and elaborated 
analytical analysis is possible2£. Here we summarize the 
most important results of this analysis for two limiting 
cases. 

In the large-size regime, L ^ Ajh or B^ 'C Bl, the 
smooth alternating deformation v(jj) has solution in the 



form of two isolated surface solitons^°. For example, near 
the edge y = such soliton solution decaying from the 
surface into the balk is given by the well-known formula 
for the sine- Gordon kink 



(82) 



tan V — tan vq exp I — 2 v2y//i j , 



where the boundary value vq can be found from the 
boundary condition ([5D|) leading to tan (2uo) — V2 cosa. 
Using this solution, one can find the surface energy and 
surface current for the edge y = as functions of the 
lattice displacement a. 



is(") 



1 
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(1 - V2 + cos2a) , (83) 

^^"2" (84) 



V2h V2 + cos2a' 



The 2a periodicity of these results is a consequence of 
the triangular lattice structure: the change of a by tt 
corresponds to the vertical lattice displacement by one 
layer. Similar solution is realized at the opposite edge 
y — L. Its energy and current can be obtained from the 
above results using the substitution a ^ a + hL. For a 
wide stack one can neglect the interaction between the 
solitons and the total Josephson current is given by the 
sum of two independent surface currents, 

J (a) ^ is{a) + js{a + hL). 

The critical current J^ can be found as maximum of J{a) 
with respect to a giving the following result in real units 



Jc{B) = J J 



$n 



2TidLB. 



-T 



2TrdLBx 



(85) 



where Jj — jjLL^ is the maximum Josephson current 
through the sample at zero field, and oscillating func- 
tion T{x) has period tt and in the range < x < 7r/2 
can be approximated by the following formula, ~^(x) ~ 
0.128 + 0.888 cos(x) + 0.021 cos(3x) ■ We can see that, 
in this regime the product B^Jc has periodicity of half 
flux-quantum per junction and reaches maxima at points 
$= dLB, =j$o/2 with S,Je,„iax « 1.035 Jj$o/(27rdL). 
This corresponds to the low-field part of the plot in Fig. 
[TBI All other properties of the sample should also oscil- 
late with the period of half flux quantum. Such oscilla- 
tions of the flux-flow resistivity in BSCCO micro-mesas 
have been first detected experimentally in Ref. \4^ and 
later confirmed by several experimental groups. 

In the small-size regime L < Ajh or B^ > Bl the 
interaction with edges dominates. As a consequence, ex- 
tended regions of the rectangular lattice appear in the 
phase diagram, see Fig. [151 The energy of the rectan- 
gular lattice, V = ±7i'/4, coincides with the well-known 
result for a single junction 
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(86) 
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and has minimum /roct — — 2 |sin(ft,L/2)| //i at a = 
—hL/2 + 6Tr/2 with 6 = sign [sin (/iL/2)]. An accurate 
analysia^S shows that the rectangular lattice is stable 
with respect to small deformations at a = —hL/2 + 7r/2 
in the regions \hL/2iT - (fc + 1/2)| < 1/4 only if the in- 
equality 



Isin (hL/2\\ < tan 



(V2L/h 



(87) 



is satisfied. These regions are plotted in the phase di- 
agram, Fig. [TS] This means that the rectangular lat- 
tices first appear in the ground state at points hL = 
{k + l/2)27r for L/h < h ^ arctan (^2) /V2 w 0.675. 
This corresponds to the dashed line shown in the phase 
diagram of Fig. [151 If, however, L/h is only slightly 
smaller than this value, the rectangular lattice becomes 
unstable with increasing current and the configuration at 
the critical current still corresponds to the deformed lat- 
tice. The accurate analysis shows that there is another 
typical value of the ratio L/h, L/h = I2 ~ 0.484, be- 
low which the rectangular lattice remains stable up to the 
critical current. 

In the region h ^ L the rectangular lattice is realized 
in the most part of the phase diagram except narrow 
regions in the vicinity of the integer-flux quanta lines, 
hL/2n = $/<i>o = k where the interaction with edges 
vanishes. Switching between the rectangular and trian- 
gular lattices in the ground state occurs via a first-order 
phase transition^*^ at the transition fields which are de- 
termined by equation 
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At high fields the critical current approaches the clas- 
sical Fraunhofer dependence for a single small junction, 
Jf{^) = >/j| sin(7r$/$o)|/k<I'/$o|- Two important de- 
viations persist at all fields and sizes: (i) Near the points 
$ — fc$Q, due the phase transitions to the triangular lat- 
tice, the critical current never drops to zero and actually 
always has small local maxima; (ii) Away from the points 
$ — fc$o the critical current is reached at the instability 
point of the rectangular vortex lattice and it is always 
somewhat smaller than the "Fraunhofer" value Jf(<&). 

In the region B ~ B^ the crossover between the two 
described regimes takes place. In the oscillations of the 
critical current this crossover manifests itself by breaking 
the $0/2 periodicity, the maxima at the half-integer flux- 
quantum points ^ = (k + l/2)$o progressively become 
larger while the maxima at the integer fiux-quantum 
points <& = fc<l>o become smaller. This crossover behav- 
ior of the critical current is illustrated in Fig. [Tni Such 
behavior was indeed observed experimentally in very nar- 
row BSCCO mesasi^'i^'iS. 



VI. THERMAL FLUCTUATIONS 

In this section we consider thermal fiuctuations effects 
for the Josephson vortex lattice. Confinement of the vor- 



tex cores in between the layers leads to strong suppres- 
sion of the vortex motion across the layers, which can 
only occur via formation of kinks. Therefore, as a first 
step, one can neglect these energy-costly displacements 
and consider only planar fiuctuations of vortices along 
the layers. This simple model describes fiuctuation be- 
havior in the most part of the field-temperature phase 
diagram but it occurs to be insufficient for description of 
the melting transition of the lattice. In general, thermal 
effects for Josephson vortices are much weaker then for 
pancake- vortex lattice and phase transformations are ex- 
pected only in the vicinity of the transition temperature. 
On the other hand, due to the intrinsic pinning poten- 
tial and involvement of the kink excitations, the overall 
behavior near the melting line is rather complicated and, 
in spite of quite extensive theoretical effort^®'^^'^'^"^^ and 
numerical simulation o 1 , there is no clear consensus 
on the nature of melting transition and structure of the 
phase diagram for magnetic field aligned with direction 
of the layers. 



A. Thermal effects for dilute Josephson vortex 
lattice: problem of intermediate phase 

A standard first step to study thermal fiuctuation ef- 
fects is to evaluate the mean-squared local fiuctuation 
displacement from the elastic energy (15^1) .^ 
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(2^)3cii(k)fc2+C44(k)fc2+c66fc2 



(89) 



Introducing reduced wave vector k as 
kx = kBzkx/y/l, ky = kBzky/y/j, k^ = kBZy/jk^, (90) 



where fcez — \/4:TtBx/^o is the average wave vector of 
the Brillouin zone, we rewrite this integral in a more 
explicit form. 



^~-- e, 




with j3 ^ 1. Evaluation of this integral leads to the fol- 
lowing result 



0.12T 



al boy^ln{bo/d)eo' 



(91) 



where qq — \/y^o/Bx and bo — a/^o/tB^ are the typi- 
cal lattice constant in the y and z directions. From this 
result one can obtain estimate for the typical tempera- 
ture at which fiuctuations become strong?^ 



Tf ^ boVl^{bo/d)eo{Tf). 



(92) 
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Unfortunately, this temperature is located very close to 
Tc where one can not use approximations behind Eq. 
([59| . e.g., neglect thermal activation of kinks and an- 
tikinks. We can conclude that the planar-fluctuations 
model given by the elastic energy (j59p is not sufficient 
to describe the melting of the Josephson-vortex latticed 
The temperature scale (|92p is much higher than the 
corresponding temperature scale for the pancake vortex 
latticei^ meaning that thermal-fluctuation effects for the 
Josephson vortex lattice are much weaker than for the 
pancake vortex lattice. 

We can estimate the typical temperature above which 
kink formation strongly influences the fluctuation dis- 
placements of the vortex lines. In the isolated line the 
typical distance between thermally excited kinks is given 

by 



crystal 



smectic 



-^kink = Ckinkexp(£'kink/r), 



(93) 



where i^kink ~ rf^o ln(7c'/Cab) is the kink energy. Usually, 
it is assumed that the preexponential factor ^kink is of the 
order of the in-plane coherence length ^ab-^^- Analysis of 
fluctuations of the order parameter near the core-- gives 
somewhat more accurate estimate ^kink ^ S,a.h\/T/deQ. 
Typical k^ contributing to the fluctuation displacement 
(|M| can be estimated as kr^ ~ Tr/fog. Therefore, the kinks 
start to contribute to thermal wandering if Lkink < &o- 
This gives estimate for typical temperature 



Tkink — -Ekink/ln(6o/Ckink)- 



In the limit 7 > ^ab/rf, we obtain 



Tkink — C?eo(Tkink) 



ln(6o/6ink) 



(94) 



(95) 



One can see that even though this temperature is smaller 
then Tf (|92|l . it is also located close to the fluctuation 
region near Tc and very slowly decreases with increasing 
magnetic field. 

The planar-fluctuations model belongs to universality 
class of the three-dimensional XY model meaning that 
the phase transition described by this model has to be 
continuous. In spite of insufficiency of this model, this 
suggests that the melting transition for the magnetic 
field applied along the layers may become continuous for 
sufficiently high anisotropy. It was indeed observed ex- 
perimentally by Kwok et al^ and by Gordeev et al^ 
that the melting transition in YBCO becomes continuous 
when magnetic field is aligned with the layers. Continu- 
ous melting of the Josephson vortex lattice also has been 
observed in numerical simulations by Hu and Tachiki2S. 
The simulation parameters in this work, however, corre- 
spond to the regime of dense lattice, which will be con- 
sidered below. 

A description of the fluctuating Josephson vortices 
taking into account kink-antikink formation is much 
more complicated problem and possibilities for analyt- 
ical progress are quite limited. General scenarios of 



liquid 



FIG. 17. Possible phases for field applied along the layers. 
Grey level illustrates average vortex density. In the interme- 
diate smectic phase suggested in Ref. [TQ density is modulated 
only in the direction perpendicular to the layers. 



Josephson- vortex-lattice melting have been discussed by 
Balents and NelsonS. They argued that an aligned lat- 
tice may melt via an intermediate smectic phase, in which 
the average vortex density is modulated only in the di- 
rection perpendicular to the layers but no order is pre- 
served in the direction of the layers, as illustrated in 
Fig. [T71 The period of density modulation has to be 
equal to the integer number of layers. The developed 
Landau theory of the liquid-to-smectic transition sug- 
gests that this transition has to be of the second order. 
Static and dynamic properties of the intermediate smec- 
tic phase have been described in detail. In particular, 
Balents and Nelson argued that this phase is character- 
ized by finite but very large tilt modulus, corresponding 
to very small transversal susceptibility ^z = B^/Hz, and 
by very small in-plane resistivity. Both these properties 
appear due to the thermally-activated "superkink" exci- 
tations, in which one vortex is moved across the layers 
by one smectic period. While the density modulation re- 
mains static and oriented parallel to the layers, these ex- 
citations may facilitate tilting of the magnetic induction 
with respect to the layers and flux motion in the z-axis 
direction. In spite of its physical appeal, the theory of 
Balents and Nelson is not quantitative. It does not pre- 
dict locations of the transitions in the field-temperature 
plane, their thermodynamic signatures, and the width 
of the intermediate-phase region. The very existence of 
the intermediate smectic phase has been not rigorously 
proven. Alternatively, the crystal may melt directly into 
the liquid via a first-order phase transition. 

More quantitative study based on the density- 
functional theory has been performed recently by Hu, 
Luo, and Mai^l. The intrinsic pinning potential in this 
study has been modeled by the cosine function and its 
strength has been used as an adjusting parameter. It 
was found that the smectic phase exists for sufficiently 
strong periodic potential only for one type of aligned lat- 
tice, which in our notations corresponds to (jn, n) — (1, 0) 
and with one empty layer between the layers filled with 
Josephson vortices, i.e., with N — 2. According to anal- 
ysis of Sec. IIVBI such lattice is realized in ground state 
within field interval [0.8 - 0.98]$o/(27r7d2). The melt- 
ing scenario via the intermediate smectic phase is most 
probable in this field range. 
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B. Elimination of the lattice rotational degeneracy 
by thermal fluctuations 



energy ([48]) , we represent the orientation-dependent part 
of the total free energy at finite temperature in the form 



The dilute lattice at small field is approximately degen- 
erate with respect to elliptic rotations, as it was discussed 
in Sec. IIVI This degeneracy is partially eliminated by in- 
trinsic pinning potential and by the corrections to the 
intervortex interactions due discreteness of the layered 
structure. The latter effect becomes noticeable only at 
high magnetic field approaching the crossover field. As 
the Josephson vortices mainly fluctuate along the layer 
direction, the fluctuation correction to the free energy de- 
pends on the lattice orientation with respect to the layers 
and also eliminates the elliptic degeneracy. Therefore the 
Josephson vortex lattice at small fields gives a physical re- 
alization of a system in which the ground state is highly 
degenerate at zero temperature and this degeneracy is 
eliminated by thermal fiuctuations. Similar behavior is 
realized in some frustrated magnetics and is known as 
"order as an effect of disorder"—. As a natural way to 
prepare the ground state is to cool system in fixed field, 
it is important to understand how the ground-state con- 
figuration evolves with the temperature. 

In this section we consider the orientation-dependent 
entropy correction to the free energy. This will allow 
us to trace evolution of the ground-state configurations 
with increasing field at finite temperature. Qualitatively, 
fiuctuations favor soft lattices, with smaller elastic con- 
stants. One can expect then that the entropy correction 
favors the aligned lattice (1,0), because for this lattice 
the shear deformations take place along the closed-packed 
direction. 

The orientation-dependent entropy correction is deter- 
mined by the short- wavelength lattice deformations and 
the long- wavelength elastic approximation of the previ- 
ous section is not sufficient. The elastic energy for planar 
deformations in the whole Brillouin zone is given by 
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where Q — {Qy,Qz) are the reciprocal-lattice vectors. 
The fluctuation correction to the free energy is given by 
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(98) 



Calculation of this correction is described in detail in 
Appendix [C] Combining result of this calculation with 
the London-limit presentation of the lattice interaction 



Sf^ = ^-[GL-^\r-^sA. (99) 

$0 7 \ £o V 7r$o J 

Numerically computed orientation-dependent correction 
ga{d, h) in the range 0.001 < /i < 0.1 is well described by 



gai0,h)^ge{h)cos{69) with ge{h)' 
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The fluctuations give the largest negative contribution 
for 9 = 0, meaning that they indeed favor the aligned 
lattice (1,0). 

Let us compare the orientation-dependent entropy cor- 
rection with the correction due to the layered structure 
considered in Sec. IIVBI We can see that these corrections 
compete: the first one favors the (1,0) orientation while 
the second one favors the (1,1) orientation. The entropy 
correction decays with decreasing fields as ^/B^ and al- 
ways exceeds at small fields the "layeredness" correction, 
which decays as B"^. We estimate that the "layeredness" 
correction exceeds the fiuctuation correction when B^ ex- 
ceeds the temperature-dependent field scale 
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with C't ~ 2.6 • 10". 



C. Fluctuations and melting of dense Josephson 
vortex lattice 

Using elastic energy (|75|) we can evaluate a mean- 
squared fluctuation of the in-plane phase 
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with Sz{kz) = sm(kzd/2) and the lattice displacement 
Un == -A.7 (i;„+i - Vn) /h 
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Renormalization of the effective coupling is determined 
by the following average 
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All above integrals are logarithmically diverging at large 
fc||. This divergency has to be cut off at fcy ^ l/^ab- As 
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usual for quasi-two-dimensional systems, the weak inter- 
layer coupling cuts off logarithmic divergency at small fcy . 
Evaluating the integrals, we obtain 

{€) « 2;^ In (Aj VCab) ; {u') « ;^ In (A,7 VCab) ; 

((x'n-i + x'n+i - 2vr,f) ~ 6 (cj^l) « ^ In (Ajh/U) ■ 

Fluctuations become strong and harmonic approximation 
breaks down when ( (w„_i+w„+i — 2u„) \ ^ 1, corre- 
sponding to {(pn) ~ 1/6 and (w^) ~ a^/3 with a = Aj/h 
being the in-plane lattice constant. This gives the tem- 
perature scale 
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As ^0(0) > T; (typically for BSCCO £^0(0) ^ 250 - 
300i^), this temperature scale usually corresponds to 
temperatures close to Tc. It is somewhat lower than 
the corresponding temperature scale for the dilute lat- 
tice (|M|) and even smaller than the temperature scale for 
kink formation in the dilute lattice (|95l) . 

We turn now to discussion of the melting transition 
of the dense lattice based on the energy ([7^ describing 
weakly coupled two-dimensional systems. Behavior of 
such system has to be similar to layered XY model®^ and 
to a layered superconductor in zero magnetic field'^^. In 
the ordered phase of such systems, below the Kosterlitz- 
Thouless temperature for a single layer, a weak inter- 
layer coupling is always relevant, can not be treated as 
a small perturbation, and it restores three-dimensional 
long-range order. The transition in such system is ex- 
pected to be continuous and to occur slightly above the 
Berezinskii-Kosterlitz-Thouless transition of an isolated 
layer given by equation T^t — '^Eq(Txt)/'2- This is in 
spite of the fact that the interplane fluctuations actually 
become strong at the temperature (jlOip which is signif- 
icantly smaller than the transition temperature in the 
isolated layer Tkt- 

Hu and Tachiki studied numerically melting transition 
of the dense lattice using the frustrated XY model^^. 
They claimed that the melting transition is continuous at 
high field and it changes to a first-order transitions when 
the field drops below B = $0/2^37^^ w 1.8<^o/'2TT-/d'^. 
It is not clear how universal is this field. In principle, it 
may be sensitive to the kink energy, which depends on 
the ratio 7d/Cab- 

Experimentally, indication of the melting transition in 
the dense-lattice regime was found in small-size BSCCO 
mesas by Latyshev et al^ exploring the temperature de- 
pendence of magnetic oscillations discussed in Sec. IV Dl 
It was found that in the field range 0.6-0.8 tesla the 
magnetic oscillations of the fiux-flow voltage rapidly de- 
crease with increasing temperature and are completely 
suppressed by thermal fluctuations at temperatures ~ 4K 
below the transition temperature. 



VII. SUMMARY 

In this review we considered in detail the static prop- 
erties of the Josephson vortex lattice following from 
the Lawrence-Doniach model in London approximation, 
which mostly describes properties of superconductors in 
terms of the distribution of the order-parameter phase. 
We reviewed the properties of an isolated vortex as well as 
the structure and energetics of the vortex lattice both in 
dilute and dense regimes. In addition to standard prop- 
erties, our consideration includes quite subtle nontrivial 
effects, such as the influence of thermal fluctuations on 
the orientation of the vortex lattice. Note that we did 
not touch on dynamic properties of the lattice which has 
became a separate large field. 
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Appendix A: Calculation of the nonlocal line-tension 
energy of a single line 

For deformations with wave vectors \kx\ 3> 1/Ac 
screening effects can be neglected and the energy vari- 
ation is determined by the phase part of energy, which 
we write using scaled in-plane coordinates, {x^y) — 
{x/jd,y/jd), 
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where 0„ (y) is the straight-vortex solution. The phase 
of the deformed vortex obeys the following equation 



t>n+l - 4>n) -Sin( 



= (A2) 



with the condition (j)i{x, u{x)) — (j>o{x, u{x)) = tt defining 
the vortex core and u{x) = u{x)/jd. In the elastic limit, 
\du/dx\ <^ 1, at distances, smaller that the typical wave- 
length of deformation, the phase approximately can be 
represented as 

0„(x,y)«0i°)[y-z2(S)]. 

On the other hand, at large distances we can use Lon- 
don approximation of Eq. (|A2p and find the phase using 
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Fourier transform. This gives for the phase perturbation 

0(1) (k) =(/)(k)-0(O)(k) 
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with {kx,ky,kz) = {ldkx,^dky,dkz) and fc^ = A;^ + 
A:^ + A:^. We wih use this result in mixed {kj;,y,z)- 
representation, which is obtained by the reverse Fourier 
transform of the above equations with respect to y and z 
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with f = (y, z). 

We spUt the total energy loss given by Eq. (jAl|) into 
the a:;-gradient and transverse parts, SF ~ F^ + F^y. The 
x-gradient part, 

■^^ "" ^ X! / ^* / ^^ (Vs(/>„)^ , 

can be computed by introducing intermediate scale, 1 ^ 
i? <C 1/kx, which splits the integral into the two contri- 
butions, from small and at large distances. The contri- 
bution from f — \J\p- -|- z^ < R, with z = n— 1/2 is given 
by 



In the transverse part, 
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we replace (pn {y, z) with (pn {y — u{x), z) and represent 
4>n{x,y) as 4>n{x,y) = 4>n\y - u{x)) + (i)n{x,y) where 
the Fourier transform of (j^nixjy) at small wave vectors 
is given by 
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We will see that the main contribution to Fxy comes from 
the distances of the order of typical wavelength of defor- 
mations far away from the core. Therefore we can expand 
with respect to </>„ and can use linear and continuous ap- 
proximation 
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with y„ = ^i?2 _ („ _ 2/2)2. The quantity 
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is determined by exact phase distribution in the core. 
Using the accurate numerical solution, we estimate Cy ~ 
0.93. The contribution from the region r > i? is com- 
puted using Eq. (|A4|) . 
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with 7e ~ 0.5772 being the Euler constant, we obtain 
Combining the parts -Flc,< and Fx^^, we obtain 
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Substituting (/"(k) and computing integral with respect 
to ky and fc^, which converges at ky, kz ^ kx, we obtain 
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Finally, combining (jASp and (jA6p . we obtain the line- 
tension energy of the Josephson vortex (j3ip which is pre- 
sented already in the real coordinates and the numerical 
constant is given by Ct = 2 exp(— 7e + Cy). 



Appendix B: Discrete and nonlinear corrections to 

the Josephson vortex phase and energy at large 

distances from the core 



The phase distribution in the JV core (j)n{y) obeys 
equation (l20l) . We will measure the in-plane coordinate 
y in units of the Josephson length Aj = jd defining the 
dimensionless coordinate y = y/jd and rewrite (PH)) in 
the form 
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4-sin [(/)„+! (y) - (/)„ (2/)]+sin [(/)„_i {y) - 0„ (y)] = 0. 



At large distances from the core, ti^ + y'^ ^ l,this 
equation transforms into the isotropic London equation 
V20 = 0. In this region 0„ (y) can be approximated by 
a continuous function (y, z) with n — >■ z. Using Taylor 
series for the difference (j){y,z + 1) — (f) (y, z), we obtain, 



sin [0 (y, z-)-l)-0 (y, z)]-|-sin [0 (y, z-l)-0 (y, z)] 

"* 5i2 "*'T2 9i4 ~ 2 laij 9i2"*'"" 
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Therefore, the phase equation to 4-th order in the gradi- 
ent (which is small at large distances) is given by 



Qy2 Qg2 



1 d^(t) 
12 9i4 



dz 






0. 



(Bl) 



This equation can be solved iteratively. For the Joseph- 
son vortex located at y = in between the layers 
and 1 the zero-order solution 0" (correct to second or- 
der in the gradients), is given by the angle (fP{y,z) — 



-tan ^[{z — 1/2) /y] (note that for z 
0O(zj/7rf,n) = c/.J-(y)in(l2lD). 
S(j)^ {y, z) obeys the equation 



n we have 
The first-order correction 



1 d*(f>" 
12 dz* 



dz 



^2j,0 



9z2 



_ 2 sin(2(/.0)-|-5 sin(4 0O) 

where f^ = y^ + [z — 1/2)^. Using the solutions of the 
inhomogeneous Laplace equations 



V^(^ = 



sin 20° 



Ini 



sin 4 



-> 



we build the solution for 
correction 

sin(2 0") 



[r, 



sin2(/)° 
4f2 

12r2 
fp) and arrive at the 



'(y,z) 



16 r^ 



5sin(4 0°) 

96 7=2 



(Inf + C54,) + ^^7;:," ' + 0(l/f4). 

(B2) 

Here we have added the solution sin(2 0°)/r2 of the ho- 
mogeneous Laplace equation with an unknown numerical 
constant C^^. Comparison of these asymptotics with the 
full numerical solution gives Cscji ~ 4.362. The result 
(IB2p is given in unsealed coordinates in (P5|) . 

In a similar way, we can derive a nonlinear/discrete 
correction to the energy far away from the core. The re- 
duced energy contribution to the Josephson vortex from 
the region f < X^b/d is given by 



£jv ^ dy 



dcfhi 

dy 



1 



cos 



l^n+1 



4>n) 



In the region f 3> 1 we can again apply expansion with 
respect to small gradient along the z axis which leads us 
to the following result 

2 , / ri J \ 2 



;jv ~ / d^v 

l«:f<SAab/(i 

' 24 I dzP- 



1 
2 

1 
24 



d0 
dy 



dz 



dz 



In the lowest order with respect to small gradients, this 
gives us the correction to the energy due to layered struc- 
ture 

2 /o,n\4n 



5e 



Jv = 



1 

'24 



d^ 



l<Sr<SA„b/ 



92^0 

dz^ 



dz 



(B3) 



In the case of a single Josephson vortex this formula is not 
very useful because the integral is formally diverging at 
small distances and is determined by the small-distance 
cut off. In the case of finite vortex density, however, 
generalization of this equation will allow us to obtain a 
nontrivial correction to the vortex-lattice energy. 

In the vortex-lattice case at finite in-plane field follow- 
ing the same reasoning we obtain the correction to the 
reduced energy per unit cell (j56p 



5u = 



d^Y 



1 
'24 



az2 j ~24 



dz 



-hy 



(B4) 



where integration is performed over the unit cell and 
4>^{v) is the vortex-lattice phase within London approx- 
imation. To estimate the dominating contribution, we 
use the circular-cell approximation for the lattice phase. 
In this approximation supercurrents flow radially within 
the cell r < ac = \/2//i and vanish at its boundary, so 
that the gauge-invariant phase gradient is given by 

- -7i — = - - + ^T for < r < flc 
r da r a^ 

where a — tan^^(z/y) is the polar angle meaning that 
dcffi/dz + hy ^ — cos(ck)(l/r — r/al). The integral is 
formally diverging at small distances. This divergency, 
however, is due to the vortex-core energy. To find the 
nontrivial correction to the lattice energy we subtract the 
diverging single-vortex term. The dominating contribu- 
tion to the rest part is coming from the second [nonlinear] 
term 



6u ; 



1 

"24^ 



da / rdr cos (a) 



1 

f a 



f Y 1 



and calculation gives the following result 



Su 



32 ^ h ' 



(B5) 
(B6) 



From fit of the numerically computed energy to this for- 
mula we obtain numerical constant Ch ~ 110. Note that 
this correction does not depend on lattice orientation 
with respect to the layers. Interaction with the layers 
also eliminates the "elliptic" rotation degeneracy of the 
lattice described in Section llVl The expansion (jB5|) . how- 
ever, is not sufficient to find the orientation-dependent 
correction to the energy. To obtain such correction one 
has to obtain the next-order expansion with respect to 
the gradients (6-th order terms). 



Appendix C: Calculation of the 

orientation-dependent fluctuation correction to the 

free energy 

In this appendix we present calculation of the entropy 
correction to free energy (|55|) based on the planar elastic 
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(a): 



_^^b: 



S 







K 



FIG. 18. ie/f pari; Josephson vortex lattice in reduced coor- 
dinates rotated at finite angle 6 with respect to the layers in 
real space so that the layers align with the crystallographic 
direction (3,1). Right part: The corresponding reciprocal lat- 
tice and illustration of two selections for the basis used in 
the calculation of the entropy correction: the basic wave vec- 
tors aligned with the lattice, Gi,2, and the basic wave vectors 
aligned with the layers, Li,2. 



energy ([M|) . To facilitate calculations we again introduce 
the reduced wave vectors, k, defined in Eq. ([90l) and the 
corresponding reciprocal-lattice vectors Q = iQy,Qz)- 
In this presentation the reciprocal lattice becomes a regu- 



lar triangular lattice with the unit vector Qq = W27r/v3 
and area of the Bravais cell equals to tt. Using new vari- 
ables, we represent $jvL(k) in the compact reduced form 
as 



*JVL(k): 



47rA2 



/.jVL(k)=^ 

Q_ 



0JVL(k) 



f^y tcjfi 



6x'- 



Q 



b^^ + Q^ 



(CI) 

(C2) 



where bx = 47rAabAci?K/$o = 2 (Xah/d) h ^ 1 and 

We will assume that the lattice is rotated at a finite 
angle, 9, with respect to the layers selected in such a 
way that the layers are aligned with one of the crystal- 
lographic directions as sketched in Fig. [181 This means 
that the lattice, in general, has the form of a misaligned 
lattice sketched in Fig. |B^ and it is characterized by the 
aspect ratio, r = ^b/a^ and the shift parameter, q. For 
computing the sum over the reciprocal-lattice vectors, 
we will use two equivalent parametrizations illustrated 
in Fig. [THl The first parametrization uses expansion over 
the two basic vector of tilted lattice, Q — nGi + 7T1G2 
with m,n = 0,±1,±2.... For such expansion we can 
simply represent the component of Q along the two main 
directions of the tilted lattice, (ki, k2), shown in Fig.lTSl 



Qi 



V3 



mQo, Q2 



m 
"2" 



(C3) 



This gives Q^ = {v? + nm + w?) Q^. The {y,z)- 
components of the wave vectors are related to the (1, 2)- 
components by axis rotation. For example, for the 



component, ky, in Eq. (|C2p we have ky = cosOki + 
sm6k2. This parametrization allows us to trace natu- 
rally the dependence on the rotation angle 9. The second 
parametrization utilizes the basic wave vectors aligned 
with the layers. 



Q = TiLi + mL2 



Lr=(0,,/^ ;L2= ViF,-«J^ 



(C4) 



This basis allows us to trace easily the dependence on 
the lattice-structure parameters r and q. It also allows 
us to reduce i5iijvL(k) to a simpler form. Substituting 
presentation (|C4p into Eq. (IC2[) and taking sum over n, 
we obtain 



>.WL 



{k)=v^ J2 



sinh 



2y^TrrK{ky^my^Trr, k^) 



irrm 



K{ky-m^/Trr,kx) cosh 
sinh [2y^TrrK{my^Trr, 0)] 



2y^TrrK{ky — my^Trr, k^) 



K{m^/Trr, 0) cosh [2y^7rrK(my^7rr, 0)] — cos{2'Kqm) 



— cos 



27r(grm + fc^y^) 



with K{ky, kx) = \ bx + ky + fc^. This formula contains 
only one summation which makes it convenient for nu- 
merical evaluations. On the other hand, the dependence 
on the rotation angle here is not obvious and is hidden 
in the dependence on the parameters r and q. 

The sums over the reciprocal-lattice vectors in Eqs. 
Wf\ and (|C2p formally diverge logarithmically at large 
Q (Q). Correspondingly, the sum over to in Eq. (IC5[) 



I 

also logarithmically diverges. This divergency is due to 
the single-vortex tilt energy and it has to be cut at the 
core size, Qy ~ ^/jd. This energy has been considered 
in details in Sec. IIII Al We will split the reduced elas- 
tic matrix, 0jvL(k), mto the single- vortex, 4'svikx), and 
interaction, 0i(k), terms. 



</'JVL(k) = (f)sv{kx) + 0i(k). 
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The single- vortex term (psvikx) can be obtained from Eq. 
(|C2p by replacing summation over Q with integration 



(AnXl/Bl) <J>3v(fc.) 



ysvyKx) 



(Pq, 



Qy + ^x 



Ql 



bx^ + Q2 + A:2 5^1 + Q2 



^ ,~ , fc2 4 Q9 

0sv(fc:r) ~ — In-y-. 



for k^ <^ A/h. In the interaction term, (?!)i(k) 



(C5) 



Using Eq. ([31]) . we obtain the line-tension term in real 
units, ^sv{kx) = TT{Bx/^o)e3kl\n{Ct/jdkx) with ej = 
Eo/jd and Ct w 2.86. This corresponds to the follow- 
ing result for the reduced line-tension term, 0sv(fcx) — 

I 



6jVL(k)- 

4>sv{kx)^ the logarithmic divergency is compensated and 
the sum over Q converges roughly at Q ^ 1. In particu- 
lar, using presentation given by Eq. (ICSp . the interaction 
term can be represented as a converging sum. 



6i(k 



E 



' iky — m^/Trr] +k1 



771 — — OO 



sinh 



2y/TrrK{ky — m^/Trr^ kx) 



V 



K{ky-m^/¥r,kx) cosh 
^SU kx,{m+S/2)y/Trr~k, 

sinh [2y^TrrK(my^Trr, 0)] 



S=±i 



K{m^/Trr^ 0) cosh[2Y^7rFK(TOY^7TT, ())\ — cos{2'r:qm) 



2Y^7rrK(fcj, — TTiy^Trr, fca;) —cos 2tt (qm+kzy^j 



J2 SU[0,{m+5/2)y/¥F\ 



<5=±1 



with 



U [Kxi l^yl 



ky^bx^ + kl + kl + (-&, 1 + kl) In ( fcj, + ^bx ^ +kl + k 



Here the terms with U [...,.. ] originate from the single- 
vortex contribution 4>sv{kx) which is properly decom- 
posed to compensate the summation divergency. In spite 
of its scary look, this formula is the most suitable one for 
numerical calculations. 

From Eq. ([M]) we obtain the entropy correction to the 
free energy in reduced form 



5fT = - 



( 



T /AnB, 



2/7 V *o 



3/2 



2n 



^ In ■ 



C 



BZ 



(2^) 



ijVL(k) 

where J„„ . . . notates the integral over the Brillouin zone 
and C is dimensionless constant. Integral over kx is for- 
mally diverging. This divergency is due to short wave- 
length excitations in the vortex cores and it does not con- 
tribute to the angular-dependent correction. To separate 
the regular anisotropic correction, we subtract from the 
total free energy the isotropic single- vortex contribution 
and represent the resulting anisotropic correction as 



SfrA^) 



T 



Bx 

$0 



3/2 



9a 



with 



9a 



d^K 



BZ 



dkx In ■ 



^svyf^xj 



0sv(fca;) +'/'i(k) 



(C7) 



(C8) 



r 



This presentation is used in Eq. ([99)1 . 

The large logarithmic factor in (j>sv{kx) in Eq. (IC5[) 
allows us to obtain a useful approximate formula for 
ga- As 0i(k) ^ 1, integral over kx converges at kx ^ 
1/ ^/ln{l/h) <^ 1 meaning that for a log-accuracy esti- 
mate we can neglect /c^^-dependence of (/)i(k). Evaluating 
the integral over kx , we obtain 



9a 



V2 



VHA/h) 



d^K 



BZ 



^i(kyz). 



(C9) 



with A ^ 1. If we neglect small parameter b~^ in 0i(kyz) 
then the integral in this formula becomes field indepen- 
dent and the only field dependence of ga for ft, — >■ is 
given by the factor [ln(yl//i)]~^/2_ 

We numerically computed the reduced entropy cor- 
rection ga for different lattice orientations and reduced 
fields h. Example of the angular dependence of ga for 
h — 0.0067 is shown in the inset of figure [121 We 
found that in the range 0.001 < h < 0.1 the orientation- 
dependent part of ga can be well fitted by the formula 
(jlOOp . The dependence geih) is plotted in figure HM Pos- 
itive sign of geih) means that the fiuctuations give the 
largest negative contribution for 9 = 0, i.e., they indeed 
favor the aligned lattice (1,0). We also can see that the 
effect occurs to be quantitatively rather small, at least in 
the considered Gaussian-fluctuations regime. 
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"--. -0.9637+0.00302 cos(68 
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FIG. 19. Inset shows the example of the numerically com- 
puted angular dependence of the reduced entropy connection 
ga{0) defined by Eq. ((C8)) for h = 0.0067. Solid squares show 
results obtained using representation for fixed lattice param- 
eters r and q given by Eq. (|C6|) . This computation is done 
for layers oriented along the crystal directions {m, n) which 
are also shown in the plot. Open symbols are obtained using 
representation with the explicit dependence on the lattice ro- 
tation angle d using the expansion (|C3p . Dashed line is the 
fit to the formula go + gecos{6d). The main plot shows the 
field dependence of the coefficient ge and the corresponding 
fit in Eq. (fTOO)) . 
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